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We considered diffusion-driven processes on small-world networks with distance-dependent ran- 
dom links. The study of diffusion on such networks is motivated by transport on randomly folded 
polymer chains, synchronization problems in task-completion networks, and gradient driven trans- 
port on networks. Changing the parameters of the distance-dependence, we found a rich phase 
diagram, with different transient and recurrent phases in the context of random walks on networks 
(or with different smooth and rough phases in the context of noisy task-completion landscapes). We 
performed the calculations in two limiting cases: in the case of rapid-network-update, where the 
rearrangement of the random links is fast, and in the quenched case, where the link rearrangement 
is slow compared to the motion of the random walker or the surface. It has been well-established 
that in a large class of interacting systems, adding an arbitrarily small density of, possibly long- 
range, quenched random links to a regular lattice interaction topology, will give rise to mean-field 
like behavior (i.e., the small- world-like random links can be treated in a mean-field fashion). In 
some cases, however, mean- field scaling breaks down in "low-dimensional" small- world networks, 
where random links are added to a low-dimensional regular structure. Examples of such cases are 
the common diffusion and the Edwards- Wilkinson process, the main subjects of this paper. This 
break-down can be understood by treating the random links perturbatively, where the mean-field 
prediction appears as the lowest-order term of a naive perturbation expansion. The asymptotic 
analytic results are also confirmed numerically by employing exact numerical diagonalization of the 
network Laplacian. Further, we construct a finite-size scaling framework for the relevant observ- 
ables, capturing the cross-over behaviors in finite networks. This work provides a detailed account 
of the self-consistent-perturbative and renormalization approaches briefly introduced in Refs. [lj 
and 0]. 

PACS numbers: 89.75.Hc, 05.40.Fb, 05.60.Cd 68.35. Ct, 



I. INTRODUCTION 



Our environment is pervaded by large-scale complex 
systems, rooted in different backgrounds: sociology, 
decision-making, economy, transportation, ecology, reg- 
ulation and transport in biology, information transfer, 
and so on [!, 0, The major characteristic of these 
systems is that they do not have any regular structure, 
and they exist over complex networks. Characterizing 
these networks is a nontrivial task, and finding the most 
important quantities to classify the static structure of 
them has been intensely researched. Recent research on 
networks has shifted the focus from the structural (topo- 
logical) analysis to the study of processes (dynamics) on 
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these complex networks [fj|. 

One of the most common phenomena that occur in 
a large number of processes is diffusion, where the flow 
across the links, and consequently, the rate of change 
of a generalized density at a given node, is driven by 
the local density gradients between the nodes and its 
neighbors. Therefore, it is important to understand 
the properties of the diffusion and the related diffusion- 
operator (or Laplace-operator) on these networks. In our 
work, besides the results, we present a technique, based 
on impurity-averaged perturbation theory, which gives 
a novel tool to study processes on networks especially 
small- world networks. 

First introduced by Watts and Strogatz 0, Q , small- 
world networks were constructed to model social interac- 
tions Q. In the original construction, on a regular lattice, 
each node is connected to k nearest neighbors, then each 
link of a node is rewired to a randomly chosen node with 
probability p. Changing the value of p, one can inter- 
polate between a regular (p = 0) and a completely ran- 
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FIG. 1: Schematic construction of a small-world network 



dom, Erdos-Renyi, network (p = 1) [T3. llTj[ . The original 
construction of Watts and Strogatz was difficult to treat 
mathematically so a new construction was suggested [12] 
combining two networks: a regular d-dimensional one, 
and on the top of that an Erdos-Renyi network (see 
Fig. [T]) with parameter p. We refer to this construc- 
tion as the "p/ain" SW network in order to distinguish it 
from power-law SW networks, defined later. Well known 
classical models (originally defined on regular lattices) 
have been studied on these networks, such as the Ising 
model Gjl [ft [ll H the XY model [H[, diffusion 
IS ' ^^^^^^ EE E3» and synchronization 

In this paper, we consider a class of networks, called 
power-law small-world (PL-SW) networks [H, HE HI • 
In these networks, on top of a regular d-dimensional 
lattice, the random long-range links have a distance- 
dependent power-law probability distribution, where the 
probability of two nodes being connected is give by 



Probability (i and j connected) 



Mr 



= Pf(r). (1) 



Here, r = is the Euclidean distance between the two 
nodes in the original regular lattice, a is the strength of 
the distance-dependent suppression of the random long- 
range links, and p is the density of such long-range links 
(i.e., p is the average number of random links on a node.) 
Af is the normalization-factor of the distribution 
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(2) 



i=l 



where L is the linear size of the d-dimensional lattice. 

Both the structure of these networks [37| and diffusive 
dynamics on these networks [2l| have been studied by 
others, in one dimension. Our results provides a detailed 
exploration of the phase diagram in diffusion processes as 
the connection topology interpolates between the origi- 
nal "plain" SW (a — 0) and the purely short-range con- 
nected (a = oo) network. One expects that when the 
probability of the random links decays rapidly enough 
(a — > oo), the large scale behavior of the processes on 



them is not affected by these links. When the proba- 
bility of the random links decays slowly (a — > 0), the 
results of the original SW network should be recovered 
0, [Hj]. In between, there is a transition which is the 
focus of this paper. Our starting point is a calculation 
on rapidly-updated networks which turns out to be the 
first-order approximation of a naive perturbation expan- 
sion for the quenched system. For small p, this expansion 
breaks down and we apply a self-consistent calculation. 

Earlier numerical work plj for d — 1, in the context 
of the random-walk, suggested that a crossover from the 
transient to the recurrent phase occurs at around a ~ 2. 
Here, we compute the asymptotic scaling properties of 
the propagator analytically (and confirm it numerically 
in d = 1), and precisely construct the full phase diagram 
for all d. This work provides a detailed account of the 
self-consistent perturbative approach, briefly introduced 
in Refs. [1] and 0], and expanded numerical results. 

In our research we will consider different diffusion- 
driven processes, defined in the following sections, on 
these networks. Two different cases will be distinguished: 
the rapidly-updated network, when the rearrangement 
of the random links on the network takes place on a 
timescale much faster than the process itself, and the 
quenched network, when the random links can be consid- 
ered static with respect to the timescale of the process. In 
the network literature, the first case of rapidly-updated 
networks is often referred to as annealed networks, and 
we will refer to them as annealed in this manuscript also. 
We will denote the annealed average of a quantity by 
]...] aim and the quenched average by [...]. In many cases, 
in the quenched system, the effect of the random long- 
range links can be understood by a general criterion [3a ] , 
based on averaging over the long-range links in a mean- 
field (or annealed) fashion (see Section [Til ). For a few 
models, this criterion is violated, leading to non-mean- 
field behavior such as diffusion on a small-world network 
with a one-dimensional spatial structure. 

Such a power-law SW network can emerge in 
a randomly-folded polymer where, besides the one- 
dimensional chain of the monomers, two distant 
monomers can be adjacent to each other, facilitating ran- 
dom long-range links, because of the spatial conformation 
of the polymer [39j . For example, in three dimensions, 
for an ideal chain, the probability of two distant sites to 
be adjacent is proportional to r~ 3 / 2 where r is the dis- 
tance of the two segments, while for a polymer chain in 
a g ood solvent, this probability is proportional to r -197 
[401 ]. Note, that in the case of small- world networks gen- 
erated by such polymer configurations, the random links 
are correlated: there is a high chance to find an other 
long-range link around an existing one of the same length 
scale. In the case of annealed networks, these correla- 
tions are "washed away" by the constant rearrangements 
of the links but, in the quenched case, the correlations 
can change the universal, large scale, behavior of the 
processes on them. In our study, we will deal with un- 
corrected random links and leave the correlated case for 
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further research. 

Another field where such structure can emerge is that 
of synchronization problems in distributed computing 
33, 41] where synchronization is achieved by introduc- 
ing random communications between distant processors. 
Choosing a power-law SW network may be preferable as 
it lowers the cost associated with communications. It 
was also argued that "wiring-cost" considerations [42j 
for spatially embedded networks, such as cortical net- 
works [43] or on-chip logic networks 44 1, can generate 
such power-law suppressed link-length distribution. In 
general, wiring cost in SW networks can result to the 
introduction of power-law SW networks [42] ]. 



for all k-s, where 9 k is the normalized eigenvector and 
Afe is the corresponding eigenvalue of the operator and 



hi- 



(8) 



The i = index is reserved for the uniform zero-mode of 
the system: 



A = 



and 



(1,1,1,...,!). 



(9) 



(10) 



A. Diffusion processes 

One case in which the diffusion equation arises on a PL- 
SW network is a macromolecule randomly moving along 
a polymer chain jumping over adjacent segments with 
nonzero probability [U [3f| HjJ [4(| ; if the macromolecule 
motion is fast compared to rearrangements of the links, 
then the network may be considered as quenched, while 
if the macromolecule motion is slow compared tothe link 
rearrangements then the network is annealed 47}. The 
equation describing random walk processes is 

d t P T , (t) = ~J2 ( A ?<r» + ffA3* ) P T „ (t), (3) 



where P r < (t) is the probability of finding the walker at site 
r' at time t and q is the relative transition rate through 
the random links. Above we have used that the diffusion 
operator A r y on SW networks, embedded in regular al- 
dimensional networks, can be written as 



Ar,r< = A° r , 



q£[ 



rnd 
r,r' 5 



(4) 



where A is the regular <i-dimensional Laplace operator 
on the underlying lattice, A rjr / is the Laplacian on the 
random part of the network, and q is the relative strength 
of the relaxation through the random links. For exam- 
ple, in one dimension, A°- = 2<5;.j — <5i.j-i — 5i.j+i- The 
diffusion through the random links is 



The Green's function, G r /y (t), of the diffusion process 
is the solution of the diffusion equation with initial con- 
dition P T >(t = 0) = <5 r '. r ". Because the diffusion equation 
has only a first-order time derivative, formally its solu- 
tion can be written as 



P T ,(t = 0) 



(11) 



Using the spectral decomposition of the diffusion opera- 
tor one obtains 

e -(A° +g A™*)t = 5^e- A **j**><\&*| , (12) 



, = vl/^x]/^, is the projector to the 
subspace of \P . It is often useful to introduce the 
Laplace transform of the Green's function, G r ' ir »(w) — 
/ °° Gr'.r" (t)e~ u ' t dt. Performing the integration over the 
time domain, one finds 



G r ',r"M = (A° + qA rnd + Uj)~,\„. 



(13) 



In order to characterize the diffusion process on differ- 
ent networks, one of our focuses is on the scaling prop- 
erties of the expected number of returns of the random 
walker by time T. For one node at site r, it is defined as 



A 



rnd 



-A" l r , ifr^r', 



(5) 



where A 



rnd 



1 if there is a random long-range link con- 



necting sites r and r'. For PL-SW networks [see Eq. (JTJ], 



Probability (A r T n J = 1) 



Af\i 



p/(|r-r'|) 



Probability (A™j = 0) = 1 -pf(\r - r'| 



(6) 



For technical purposes, it is often useful to employ the 
spectral decomposition of the diffusion operator. In the 
orthogonal eigensystem: 



A** 



(7) 



F r (T) 



dtP r (t) 



dtG r ^ r (t). 



(14) 







Averaging over all nodes in the network one then has 

i-T 



,(*)• 



(15) 



We are interested in its long-time behavior in the thermo- 
dynamic limit. A random- walk process is called recurrent 
if in an infinite system, the expected number of returns 
is infinite as T — > 00 and transient if, on average, the 
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walker visits the origin only finite number of times. Us- 
ing the spectral decomposition of G r . r (t) from Eq. (fl"2f . 
one obtains 

F(T) = * f ^e-^X;i«Jl J * 

fc#0 

The above form is useful for exact numerical calculations 
for a given realization of the network. 

To analytically extract the leading-order asymptotic 
scaling properties for F(T), instead of the sharp cutoff, it 
is sufficient to use a kernel that suppresses the integrand 
for t T. To make our calculation easier, we chose 
J °° G> )r (i)e — */ r dt (48|. From the spectral decomposition 
of the Green's function, and repeating the same steps as 
above, now for the exponential cutoff, one finds 

HT) oc ^ £ G rA . = 1/D = j- d £ + 1- 

(17) 

In order to make the formalism similar to that of the 
other phenomena described in this section, let us intro- 
duce a modified GF, 

G r , r 'H = G VH--U- (18) 

It is basically the GF defined in the space orthogonal to 
the uniform zero mode. Note that, in the next section, 
defining the GF in this subspace will be necessary be- 
cause the coupling matrix is non-invertible. In terms of 
this modified GF 

F ( T )«7dE^r^-l/T) + -^. (19) 



B. Synchronization in task-completion networks 
and the Edwards- Wilkinson process 

In this section, we motivate the Edwards- Wilkinson 
(EW) process [4!| on networks, which can be thought of 
in terms of a synchronization paradigm in a noisy en- 
vironment. Consider a distributed task- processing net- 
work, where the processing nodes, in order to schedule 
and perform new tasks, must wait for the results delivered 
by other nodes. This task-dependency between nodes 
correspond to the links, which in general, can be asym- 
metric or directed, quenched or time-dependent. Exam- 
ples include distributed-computing systems [411 ]. man- 
ufacturing supply chains, or e-commerce-based services 
facilitated by interconnected servers 50]. Understand- 
ing the scaling behavior of the fluctuations in our model 



will help us better understand the generic features of 
delays and back-log formation in large-scale networked 
processing systems. In particular, for certain distributed 
computing scheme s (p arallel discrete-event simulations) 
[UIH, HU HH HH, [56j, it was shown [E3] that on regular 
grids, the evolution of the task-completion (or synchro- 
nization) landscape is governed by the Kardar-Parisi- 
Zhang (KPZ) [58| equation. The resulting landscape is 
"rough", corresponding to a large spread in the locally- 
completed tasks on different nodes. In general, the width 
provides a sensitive measure of the average degree of de- 
synchronization in task-completion networks, 

1 N 

™ 2 (*) = ^E(M*)-m*)) 2 - ( 2 °) 

i=i 

Here, hi(t) is the local progress on node i, h — 
{l/N)Y^iLih>i(t) is the mean progress, and N is the 
number of processing nodes in the system. When the 
network has an underlying d-dimensional structure with 
linear size L, one also has N = L d . 

In order to improve the uniformity of the progress of 
the nodes in the above distributed-computing example, 
ultimately leading to better scalable performance, one 
must suppress large fluctuation in the synchronization 
landscape. For example, one can construct a scalable 
autonomous synchronization schemes where the nodes, 
in addition to their nearest neighbor on the grid, also 
communicate with random (possibly distant) neighbors 
(with a very low frequency); the sole purpose of commu- 
nications through the random links is to keep the spread 
of the task-completion landscape under control [33|, [4JJ • 
In this case, the underlying synchronization network is 
a SW network (which can be quenched or annealed, de- 
pending on the implementation). 

While the precise local synchronization rules give rise 
to strongly non-linear effective interactions between the 
nodes, one can gain some insight by considering the lin- 
earized version of the effective equations of motion. As 
it was shown for the basic distributed-computing syn- 
chronization problem [571 ] . the dynamics, neglecting non- 
linearities, can be effectively captured by the EW process 
[49l ]. The generalization of the EW process to complex 
networks 0, 0] (focusing on PL-SW networks in this pa- 
per), is given by 

d t hi(t) =-J2 &ijhj{t) + Vi(t) , (21) 
i 

where Ay is the diffusion operator on the network and 
r]i(t) (without loss of generality) is a delta-correlated 
white noise with variance 2 : 

(rn(t)Vj(t'))='2Si j S(t-t'). (22) 

In general, the symmetric relaxational couplings between 
the nodes, Aij (An = 0), can be weighted. The general- 
ized network Laplacian then can be written as 

Ay = $ij ^ ' Aik — Aij . (23) 
k 
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For unweighted networks with unit coupling strength Ay- 
is simply the adjacency matrix. 

Clearly, the precise scaling behavior of synchronization 
landscape (its width, in particular) is strongly affected 
by both the nonlinearities and the underlying interaction 
topology. Our results provide a detailed account on the 
network's effect on the EW synchronization problem. Be- 
cause of the linear nature of the model, this is a much 
simpler problem than the actual synchronization dynam- 
ics, but, as we shall see in the following sections, the be- 
havior of it is still nontrivial, particularly on quenched 
SW networks. The understanding of the EW process (on 
networks) is a necessary first step before tackling more 
complicated, nonlinear problems. 

Using the Fokker-Planck description of the EW pro- 
cess, it can be shown [59( that the steady-state proba- 
bility distribution function (PDF) of the possible surface 
configurations is 



P[h] oc exp(-H[h 



(24) 



where h = {hi}f =1 is the configuration-vector and H[h] 
is the steady-state Hamiltonian[81| of the surface 



H[h] = i ^2 hiAijhj 



2 ^ 



Aij [hi 



(25) 



The steady-state Green's function (GF) or propagator 
of the system is defined as the two-point correlation func- 
tion. Since the steady-state PDF is Gaussian, the two- 
point correlation function is the inverse of the coupling 
matrix. In our case, A is non-invertible because of the 
uniform cigcnmode with zero eigenvalue [Eq. (JTDJ)] . 
The singularity of this mode can also be traced back to 
Eq. (I21|) , in that the mean height performs a simple ran- 
dom walk and has a diverging variance in the limit of 
t — > oo (for any system size), i.e., its steady-state vari- 
ance does not exist. Equivalently, the steady-state PDF 
Eq. ([M]) with the Hamiltonian Eq. (J23J) is ill-defined, 



</> 2 



u oo". This problem can be overcome by con- 



structing observables which lie in the space orthogonal 
to the zero mode of the network Laplacian. It can be 
achieved by, e.g., measuring height fluctuations from the 
mean [see Eq. (|20p]. Hence, the appropriately defined 
GF is 



Gij = ((hi - h){h,j - h)) = A- 



(26) 



where A -1 is the inverse of A in the space orthogonal to 
the uniform zero mode. Particularly useful for numerical 
purposes, one can employ the spectral decomposition of 
the inverse of the network Laplacian in the space orthog- 
onal to the zero mode [821 



G, 



(27) 



k=l 



Note that since the GF is constructed from the height 
fluctuations measured from the mean, the zero mode au- 
tomatically drops out from the sum. 



Employing the GF, the average width can be expressed 



2 >=(^i>-^) 2 )=^!> 



1=1 



or, for numerical purposes, 



(w 2 ) = 



L«-l 

— T — 

L d ^ A fe 

k=l K 



(29) 



C. Other examples and connections between 
transport, resistor networks, and the EW 
synchronization problem 

In the network research literature, an "abstract" trans- 
port efficiency can be considered in analo gy w ith electric 
conductance properties of a network [oOl l6ll. loH . lo3l . loll . 
[6o1 . [66| . The average resistance of a network gives a sim- 
plified but palpable interpretation of its transport prop- 
erties. If the average network resistance is low, it can be 
interpreted as a network with good transport properties; 
if the resistance is high, the network is likely to be a poor 
choice for transportation purposes. 

Consider an arbitrary network with N nodes where 
the conductance of the link between nodes i and j is 
A^. Ski pping the details (which can be found in, e.g., 
Refs. (60L 164 l66l|). using Kirchhoff's and Ohm's laws, 
one can show that the effective two-point resistance of 
a network can be expressed in terms of the same net- 
work propagator (GF) as defined in the previous exam- 
ples [Eq. (USD], 



R, 



"|- Gq 



2G, 



(30) 



Comparing the above formula to those of the EW pro- 
cess on the same network (and with the same set of 
{Ay}), connections between the two different problems 
can be readily obtained [H, [66| . Calculating the height- 
difference correlation function, 



Gjj 2Gij — Rij , 



(31) 



is equivalent to measuring the two-point resistance on the 
same network. The average resistance (averaged over all 
pairs of nodes) becomes 



R 



N(N-l)^2 RlJ ~ 



2N 



2(w z 



(32) 



which is twice the average width for sufficiently large 
networks [Hj]. These simple relationships are useful in 
understanding intuitively how the network topology in- 
fluences the behavior of different processes in networks. 
Also, note that the above relationships, Eqs. (I?!]) and 
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(|32p . are valid for an arbitrary graph. For networks em- 
bedded in a d-dimensional regular grid, such as the one 
considered in this paper, N = L d . For example, for a 
one-dimensional chain with unit resistances between each 
adjacent nodes, the two point resistance is Rij = \i — j\, 
therefore, without any calculations, the height difference 
correlation function is also ((hi — hj) 2 ) = \i — j\- 

There is also an extensive literature on the connections 
between resistor networks (with link conductances Aij) 
and discrete-time random walks (RW), with transition 
probabilities Rj = Mjl^Z^Mk from node i to node j 
[HI, [H, [69|, [z3|. Our results, to be obtained in this 
paper for the GF on distance-dependent SW networks, 
hence provide answers for these related problems. 

For example, in transport or flow problems on net- 
works, one often considers the "load" of the nodes or the 
links, which can be ca ptu red by an appropriately defined 
betweenness measure IZllzlIzi (for a given dynamics). 
This observable captures the amount of traffic, informa- 
tion, or data, the nodes or the links typically handle in 
the corresponding information, social, or infrastructure 
network. For example, in random routing or search prob- 
lems on complex networks, an appropriate node between- 
ness can defined as the expected number of visits to node 
i by a random walker (starting at a source node and be- 
fore reaching its randomly chosen target) , averaged over 
all source and target pairs 0, [75| . It can be shown to 
be given by [66| 



T,k A i 



■R, 



(33) 



where R is the average two-point resistance of the asso- 
ciated resistor network with link conductances A y - . The 
node betweenness can be used to obtain the congestion 
threshold in queue- limited network traffic [HI, [7J, [T5l ] . 
Similarly, one can show that the link betweenness bij for 
the same RW problem (the expected number of times a 
walker travels across the link between nodes i and j, av- 
eraged over all source and target pairs) can be written 
as 



bij — AijR . 



(34) 



The link betweenness, capturing the load of the connec- 
tions between the nodes, can be used to extract the con- 
gestion threshold in bandwidth-limited traffic. Further, 
one can also show that the average first-passage time for 
the above random- walk problem (av erag ed over all source 
and target pairs) can be written as [69ll7Cj] 



R . 



(35) 



As all of the above examples illustrate, the average two- 
point resistance of a graph R (or, equivalently, the aver- 
age steady-state width of the associated EW synchroniza- 
tion landscape) plays a fundamental role in the efficiency 
of prototypical transport problems. 



II. ANNEALED SMALL- WORLD NETWORKS 

In the literature of networks the disorder of links is 
considered annealed or mean-field (MF) when it can be 
approximated by its average value. Such a situation can 
emerge when the process investigated on the network 
feels only the average effect of the links due to the dy- 
namics of the system, in other words, when the timescale 
of the changes in the configuration of the links is much 
shorter than the characteristic timescale of the process 
on it. 

If O is some observable depending on a random po- 
tential V, accounting for the effect of the random links, 
then [0] arm is calculated by replacing V with its expected 
value: 



[o(v)Y 



0([V]) 



(36) 



Specifically, for the propagator of the diffusion operator, 
the focus of the present work, one has 



i 



A° + qA rnd , 
A + q[A rnd ] 



(37) 



Both in the annealed and the quenched case, the ensem- 
ble is translationally invariant, and so will be the aver- 
ages of the observables of the system. [G r \ r "] ann will 
only depend on the difference r' — r" . Let us introduce 



G ann (r) = [G r ^ +r ] c 



(38) 



where r is the spatial separation of two sites. Since 
translationally invariant operators are diagonal in Fourier 
space, it is useful to calculate the Fourier (or fc-space) 
representation of them. In Fourier space the GF can be 
written as 



G an "(k) = 



1 



k 2 + q[A rnd ](k) 



where k 2 and [A rnd ](k) are the Fourier transform of 
A and [A rnd ] respectively. The behavior of [A rnd ] for 
l/L<Cfc^;l/ais calculated in Appendix El yielding 



p if a < d 
[A rnd ](k) cx { pk a - d ifd<a<d + 2 
pk 2 if d + 2 < a , 



(39) 



with logarithmic corrections at the boundaries of these 
different a regimes. In order to keep the treatment 
of the annealed and the quenched system in parallel, 
we also call [A rTld ] the annealed self-energy, £ aim (r) = 
[A rnd ] (r) . The real-space behavior can be obtained by in- 
verse Fourier transformation where necessary and, specif- 
ically, 



G" 



w4e 



i 



L d k 2 



q [A rnd }(k) 



(40) 
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In the present and the next sections, the behavior of 
G aim (0) will be investigated. In Section [V] physically 
measurable quantities will be associated with G ann (0), 
like the surface width in the EW model and the expected 
number of returns of random- walk processes. 

For calculating scaling properties, we take the contin- 
uum limit of these quantities. In this limit, r is no longer 
constrained to the lattice sites, so it can point in any 
direction with the condition a < |r| < L, the sums re- 
placed by integrals, and Kronecker Sij by Dirac S(r — r'). 
To calculate the scaling properties of G arm (0), with re- 
spect to p, q, and L, the sum can be approximated with 
its integral- form limit: 



C. Higher dimensions, d > 2 



G orm (0) 



l/a 



dkk 



d-l 



S d 

(2n) d J 1/L k 2 + q [A rnd ](k) 



(41) 



where Sd is the surface of a <i-dimensional unit sphere. 

Mention must be made that mapping out the different 
phases of the annealed network is a relatively easy task, 
and has been investigated by others in other contexts 

[Hill]. 

The results can be summarized as follows. 



A. Phase diagram in one dimension 

The results for one dimension are sketched in Fig.[2]and 
the details of the calculations can be found in Appendix 
|B1 The relationship of G oml (0) with measurable quanti- 
ties is explained in more detail in Section [V] For a < 2, 
the integral, Eq. (|4*Tj) . converges at both ends, so is sys- 
tem size independent for large L and diverges as p — ► 0. 
There are two regimes, distinguished by their scaling be- 
havior (reflecting on the properties of the underlying ran- 
dom walk/surface) : (1) transient/smooth phase I, where 
the system behaves as if the long-range links were uni- 
formly distributed, and G ann (0) is a-independent; (2) 
transient/smooth phase II, where G ann (0) is still finite in 
the thermodynamic limit, but has an a-dependent diver- 
gence as p — > 0. For a > 2, G ann (0) has infrared diver- 
gence resulting in two system-size dependent regimes: (1) 
recurrent/rough phase I, where the system behaves as if 
there were no long-range links, and G ann (0) diverges lin- 
early with the system-size, L; (2) recurrent/rough phase 
II, with sublinear scaling with respect to the system-size. 
At the boundaries of these phases, logarithmic correc- 
tions are present. 



B. Phase diagram in two dimensions 

The behavior of the annealed two-dimensional system 
is summarized in Fig. [3J Three different phases are 
identified: a, logarithmically smooth/transient phase when 
a < 4, a doubly-logarithmic rough/recurrent phase when 
a = 4, and a logarithmically rough/recurrent phase when 
a > 4. The details of the calculations can be found in 
Appendix [B] 



In higher dimensions, the integral Eq. (|4"Tj) diverges 
for large k-s but has a well defined cutoff because of the 
lattice spacing, a. As a result G ann (0) is always finite, 
which means that the surface is always smooth and the 
random walker is always transient. 



III. QUENCHED SMALL- WORLD NETWORKS 

In this section we will develop a perturbation expan- 
sion, based on impurity-averaged perturbation theory 
77], to calculate the quenched average of the GF over 
the ensemble of random networks. The reason for in- 
vestigating such an average is the so-called self-averaging 
hypothesis, that for large enough networks the behavior 
of the GF, or in general any global quantity, will only 
deviate slightly from the behavior of its quenched aver- 
age. This assumption is due to the recognition that if a 
random system is large enough, the GF samples through 
a large amount of randomness, which is a large number 
of random links in our case. In effect, the GF will be 
the same or very similar from one realization to an other 
and for large enough system sizes it will converge to its 
quenched value. In principle, this hypothesis can always 
be checked or proven after setting up a formalism to cal- 
culate quenched averages over the randomness. To do 
so, one has to obtain the quenched variance of the GF. If 
the variance disappears in the TD limit, the hypothesis 
is proven to be valid. 

For example, the probability distribution of the ran- 
dom walker over the nodes of the network for large times 
is expected to be self-averaging. One can imagine that af- 
ter long enough time the random walker will visit a large 
number of sites, and therefore encounter a large number 
of random links. 

In this section, we will investigate the limit when the 
density of random links is small (p <C 1). Due their small 
number, we will treat the diffusion through these random 
links as perturbation of the <i-dimensional diffusion. As 
a first step, we will consider the "pedagogical" example 
of a simpler problem, where there is only a single long- 
range link in the system, which will help us to understand 
the different terms of the more difficult problem of the 
perturbation expansion. Second, we will set up a naive 
perturbation expansion and show the breakdown of it. 
Third, we will reformulate the expansion to obtain the 
leading-order behavior of the GF. For special cases, we 
will obtain higher-order corrections and investigate their 
behavior. 



A. The example of a single link 

In this section we will consider two problems that have 
similar mathematical structure as that of the problem 
of the quenched network. First, we will calculate the 
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a: 



annealed results: 



oo 



transient/smooth 

phase I 
G(0,C0=1/T) , 
G(0) ~p-" 2 

quenched results: 



transient/smooth 

phase II 
G(0,(Q = 1/T) , 

G(0) ~p*t 



recuirent/rough 
phase II 
G(0) ~L™~ 2 , 

G(0,<D=1/T)~ 

q-2 
X ot— 1 



recurrent/rough 
phase I 
G(0)-L, G(0, G)=l/T) ~ T 1 



transient/smooth 

phase I 
G(0,C0=1/T) , 
G(0) -P" 1 



transient/smooth 

phase II 
G(0,G) = 1/T) , 

G(0) ~pli 



recun"ent/rough 
phase I 

G(0)~L, G(0,Q) = 1/T) ~ T 1 



recurrent/rough phase II yt.pi-2 
G(0) -L 71 " -2 , G(0, 0) = 1/T) ~ TtfPM 
2<v(p)<3 



a: 



FIG. 2: The one-dimensional phases. 



2 
annealed results: 



oo 



logarithmically transient/smooth 
phase 

G(0,0) = 1/T)G(0) ~lln(p)l 
quenched results: 



logarithmically recurrent/rough phase 

G(0,C0=1/T)~ ln(T) , G(0) ~ ln(L) 

doubly logarithmic recun"ent/rough 
phase 

G(0,CO=1/T)~ ln(ln(T)) , G(0) ~ ln(ln(L)) 



logarithmically transient/smooth 
phase 

G(0,(0 = 1/T) ? G(0) ~ 
lln(p/lln(p)l)l 



logarithmically recurrent/rough phase 

G(0,0) = 1/T)~ ln(T) , G(0) - ln(L) 

sub-logarithmically recurrent/rough 
phase 



G(0,G) = 1/T)~ r/ (p) (T) 



o<^ P )<i 



G(0) - lr/ (p) (L) 



FIG. 3: The two-dimensional phases. (For sake of simplicity, factors of a, the microscopic cutoff, were omitted from the 
formulas and can be calculated by dimensional analysis.) 



Green's function of a diffusion process with a regular 
d-dimensional diffusion plus a "perturbation" which ac- 
counts for particles escaping from the origin with rate q. 
Second, we will consider the Green's function of a pro- 
cess that will consist of a ti-dimensional diffusion plus 
diffusion between the sites a and b with rate q. 

For the first problem, the diffusion operator describing 
the process is 

r = A + qV, (42) 



where the perturbation to the regular diffusion is Vij = 
Sio6 j. For this process, the GF can be calculated as 



A + qV (1 + qVG°)A 
= G° — G°qVG° + G°qVG°qVG° — ... , (43) 

where G° = (A )" 1 is the GF of the Laplacian on the 
regular lattice. In Fig. |4] a diagrammatic notation is in- 
troduced to represent the above series expansion. Sin- 



9 



gle lines, double lines, and crosses denote G°, G, and 
V respectively. Using the specific form of V, G can be 
calculated as 

Gij = G®j — G i0 qG j + G%qGQ qG®j 
- G i0 qG QO qG QO qG j + ... 
= G { j - G%q (l - gGpo + ^Goo^Ggo - ...) G j 

= ^ ~ G ° M ( 1 + gG8 ) ^ (44) 

To have a better grip at what this result means, let 
us calculate G for one dimension. From Appendix [D] 
Gj?. rtsL-\i-j\. Substituting this to Eqn. |@3J) and 
noticing that = {(qL) -\ +1)qL « ^(1 + -fc), 

G« « L-|i-j|_(£-N)^(i + ^)(£-|j|) 
= J + KI + lil-K-il- (45) 

Before moving on to the next problem, in the case of 
this one-dimensional lattice, let us have a closer look at 
the perturbation expansion of Eq. ([44]) in light of the 
exact behavior of the GF, Eq. (|4"5| . Comparing the first 
and the second order term of the perturbation expansion, 

G%(q)G° 0j « G Q t0 (q 2 L)G° Q] . (46) 

Even though the second order term is higher order in q, 
it diverges with the system size. As it can be seen from 
the exact solution, this divergence is not the property of 
the full GF, G, but that of the ill-defined perturbation 
expansion. In the next section, a similar problem will 
be encountered which will be resolved by reordering and 
resumming the series expansion. 

I 



B. From a naive to a self-consistent single-link 
perturbation expansion 

In the case of the random networks, a similar pertur- 
bation expansion can be done as in the previous section 
for a given realization of the random links. The method, 
Eq. (|43[) . is still valid though in this case the perturba- 
tion, V, is more complicated: it consist of a sum of terms 
like (|4"T|) corresponding to the random links in the sys- 
tem. In principle, the GF could be calculated for any 
long-range-link configuration. Though, such an expres- 
sion is not really useful for extracting the dependence 
of the GF on the parameters of the probability distri- 



For the second problem, the perturbation is, 
q y(a,b) _ y §. a § aj „ q S ib S aj + q § ib 8 bj - qS. la 5 bj . (47) 

The meaning of the first two terms is that particles disap- 
pear at node a with rate q and appear at node b. The last 
two terms are responsible for the inverse process from b 
to a. In this case we have to solve the same equation, 
Eq. (|43p . as before only with a different perturbing po- 
tential. 

To make the calculation easier, let us rearrange the 
formula of V {a ' b \ 

q V^' b) = q(S ia - S ib )(d aj - 6 bj ) = q{UW) l3 , (48) 

where we factorized y( a ' b ) as a product of two matrices 
Uik = (Sia - S lb ) and W k j = (S a j - S bj ). This way, the 
n th (n > 1) term of the series expansion of Eq. (|4"3"|) will 
become 

G°qU(WG qU) n ' 2 WG . (49) 

The three parts of this expression are (G°U)ik = 
Gl - G%, (WG )^ = G% - G%, and (WCPU)* = 
(S ak - S bk )G° kl (S la - S lb ) = G° aa - Gl b - G° ba + G° bb = 
2(G" a — G" fc ). In the last step, we used the fact that 
G° is symmetric and translationally invariant. There- 
fore, the n th term in the series expansion is G° ik (8k a — 
5 kb )q (2q(G° aa - G° ab )) n ~ 2 (S al - 5 bl )G° lr Let us compare 
it to the n th term of the previous problem: there we had 
C- k 6koQ (^Gqo)' 1-2 SoiGij. As it can be seen, the second 
problem has the same structure in its expansion as the 
previous one. The closed form of the GF is 



(50) 

I 

bution of the links. Averaging over all the realizations 
will get rid of the details due to the specifics of a given 
realization and provide us with the dependence on the 
statistical properties. 

Here, a similar diagrammatics will be used as in the 
previous section in order to keep track of the different 
terms in the expansion. For the definition of the symbols 
see Fig. [S] 

First, let us take the disorder average, represented by 
[ ], of the perturbation expansion of the GF (see Fig. [5]). 
In this expression, the higher- and higher-order terms 
have higher- and higher-order moments of the pertur- 
bation potential, A rnd . In probability theory, moments 



' — r*° — (c° — c° \ ( 1 i ir n — c° \ 

ij ~ ij { ia lb) {l + 2q(G° aa -G ab )) ( a3 bl >- 
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G G° G° qVG° G° qVG°qVG° 

— — X -|- x X 

— X x X + . . . 

FIG. 4: The simple diagrams of the example of a singe link. 




FIG. 5: The averaged Green's function in terms of the mo- 
ments and cumulants of A rnti . Single lines, double lines, and 
crosses denote G°, G, and qA rnd respectively. [ ] denote the 
disorder average of an expression with respect to the random 
potential A rnd . Crosses connected by dashed lines represent 
the cumulant of those random variables (see Fig. [6j . 

are not really suitable for any type of expansion since 
higher- and higher-order terms have larger and larger ex- 
pectational values. Cumulants are more appropriate to 
describe a random variable since they decay as their order 
tends to infinity for well behaving PDF-s. 

From a diagrammatic point of view, an important 
property of the cumulants is that the moment of an n th - 
order scattering is the sum of all the possible combina- 
tion of cumulants of the n random potentials (see Fig. [6]) . 
One might consider this as a recursive definition of the 
77,' ft -order cumulant in terms of the lower-order cumulants 
and the n -order central moment. It is also used to de- 
fine the self-energy. In the theory of disordered systems, 
the self-energy is a potential added to the unperturbed 
Hamilton operator and accounts for the average effect of 
the impurities, i.e. the random links, in the system. The 
recursive property of the cumulants is also deployed in 
self-consistent calculations where single lines are substi- 
tuted with double lines as a result of the recursive resum- 
mation of different order single-line diagrams. 



Order of cu- 
mulant 


Binary distribution 


Poisson distribu- 
tion 


1st 


P 


P 


2nd 


P-P 2 


P 


3rd 


p-3p 2 + 2p 3 


P 


4th 


p - 7p 2 + 12p 3 - 6p 4 


P 



TABLE I: Comparison of the cumulants of the binary and the 
Poisson distributions with expected value p. 



C. Cumulants of the random perturbation 
potential 

In the case of PL-SW networks, the entries of the ma- 
trix of the perturbation potential are described by binary 
PDF-s: there is either a random link between two sites or 
there is none. These PDF-s are characterized by p: the 
probability of having a link between two arbitrary sites, 
i and j, is pf(\i — j\) and the probability of having none 
is 1 — pf(\i — j\) (see Eq. (JlJ). The cumulants of this 
distribution are shown in Table 3.1. 

Observe that all the cumulants are of order p to lead- 
ing order. The cumulants do not decay as the degree of 
them increases. Some might think that for PDF-s with 
such a property an expansion in terms of the cumulants 
should not work since all terms are of the same order 
with respect of the expansion parameter, p. Later on 
in this section we will revisit this question. As for now, 
let us assume that a cumulant expansion could work in 
principle. 

The random long-range links are independently dis- 
tributed. Defining 



(o,fc) _ J 1 w ith probability pf(\a - b\) 

[ with probability 1 — pf(\a— b\), 



the contribution of each random link to the perturbation 
potential can be formulated as 



qAVf = J2qV<f' b) x( a > b \ (52) 

a<b 



where V^ a,h ^ is the single-link diffusion operator between 
FIG. 6: The second and the third moments in terms of cu- a and b introduced in the previous section (see Eq. (07])). 
mulants. The notation is the same as in Fig. Therefore, the cumulants are: 
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a<b 

[ q ^f q ^i d ] c = E^ 6) E^rV a ' 6) * M) ]c 

= E^'^E^^K^"'^) 2 ]^^ (53) 

a<6 c<d 

[gA^gA^A™^ = E9 3 ^ 6 V^ 6 V^)[(^)) 3 ] C . 



a<fc 



The n order cumulant of A will be the function of 
that of a;( a ' h ) . To make the analytics slightly easier, let us 
change the statistical properties of the random networks 
a little in a way that does not change the behavior of the 
ensemble in the p — > limit. Instead of binary random 
variables, as introduced in Eq. (|5"l"j) , Poisson distributed 
ones will be used with 



Probability( 



(54) 



Using Poisson distribution, a negligible fraction of the 
nodes will have multiple links between them, but this 
fraction is so small (0(p 2 )) that it will not affect the 
average properties of the system. From a mathematical 
point of view, though, the forms of the cumulants became 
much simpler 



[(x^r] c =pf(\a-b\) 



(55) 



as it can be seen in Table MI CI 

Let us have a closer look at the expressions in Eq. (154")) 
and interpret the mathematical formulas. For exam- 
ple, the third order cumulant consist of a sum of terms 
like V^°' Vfci' Vdm ■ Comparing it to the terms of the 
single-link problem of the previous section, one can see 
that it corresponds to the 3 rd order scattering off a single 
link between a and b. In general, the n th order cumulant 
is a sum of n th order single- link scatterings weighted with 
the expected number of those links, pf(\a — b\). 

Using the specific form of pf(\a - b\) and a - b ) [84| . 
let us calculate the exact form of the different cumulants: 



g[A™ d ] c = E(** - *»X*«i - *«) 



a<b 



Af\a-b\ a 



qp if i = j 



(56) 



operator. Note that 



&3 



rnd~\ 
ij J 



(57) 



namely, the diffusion is " conservative" . 

For the higher-order cumulants, in order to make it eas- 
ier to interpret the results, let us have Gij in between the 
random perturbation operators, as they will appear in 
the formulas of the following sections (see Section UlI Djl . 



q 2 p 





" 5ib)(Sak 


- Sbk) 


a<b 






Gkl(Sla " 


~ 6lb)(S a j 


Sbj)' 


E( 5 » - 


- S ib )(5 aj 




a<b 






2(G aa — 


Gab)q 2 p 




TV a 


-b\ a 





Af\a-b\<- 



(58) 



This is also a superdiffusion operator as in the previ- 
ous case but the weight of the matrix-entries changed to 

2 ^ G j\f\a~b\ , * q ~ ■ ® ne can eas ily generalize to the nth order 
cumulant: If there are n crosses connected by a dashed 
line, we call these single-link diagrams, and when dou- 
ble lines are running in between them, the corresponding 
matrix is 



Erx x \rx x \(2<l{Gaa — G a b)) n 1 qp . . 
{5 ia - 5ib){5 a j - 5 bj ) — — (59) 

a<b 



Af\a-b\ a 



i.e. the first cumulant is the annealed "superdiffusion" 



resulting from the aforementioned properties of multiple 
scatterings from a single link operator, V^ a ' b \ The gen- 
eralization to cases when there are different operators in 
between the cumulants of the perturbation potential is 
straightforward. 
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[=] = © — + — © — ©- 

= ©= 

© = M - x + >> -- x + x— 



©=[x]— x^=x • *: :x- x 



— x X = X x — x y y x 



FIG. 7: The introduction of the self-energy, E. The notation 
is the same as in Fig. [5] 



D. Introduction of the self-energy, self-consistent 
resummation of the diagrams 

As in quantum field theory, the diagrams of Fig. can 
be rearranged so that the self-energy is defined. The self- 
energy consist of diagrams that are irreducible, i.e. they 
cannot be cut into two pieces by cutting one of the propa- 
gators, G°, in them (see FigUJ). In terms of mathematical 
equations: 



G = G° - G°EG = (A + £)~ 



(60) 



At this point, let us compare some higher-order dia- 
grams to the first-order one in the self-energy. Let us 
consider diagrams with n crosses connected by a dashed 
line and single lines running in between them in one di- 
mension. In this case, {G° aa — G® b ) cx ?*, where r = \a — b\. 
The contribution of such a diagram to the self-energy is 



lP (2r)' 



Mr a 



(61) 



and 



5#(0) = £5#(r)«g B p 
cx q n - x V L n - x - a - 



(2r 



,71- 1 



Afr a 



-dr 



oo, 



(62) 



for some large enough n. Just like in the example of a 
single link, the diagrams diverge with the system size. 
These divergences can be avoided by including enough 
intervening scattering events of the propagator between 
two crosses In terms of diagrams, it means that a 
set of diagrams are resummed in between two crosses. If 
there is a subset of processes with intervening scatterings 
that will avoid the system-size dependent divergence of 
the higher-order terms, then resumming over all the pos- 
sible scattering events between two crosses should also 
work since it includes those processes too. In terms of 
diagrams, it means that single lines are replaced by dou- 
ble lines (see Fig. To avoid the overcounting of the 
diagrams, one must note that, expressing E in terms of 
double lines, some diagrams of E with single lines will 
become "obsolete." For example, in Fig. [8j there is no 
double line counterpart of the fourth diagram of E in 
Fig. [7J This is so because such diagrams are already in- 
cluded in the second diagram of Fig. [51 In general, to 



FIG. 8: The self-energy self-consistently. The notation is the 
same as in Fig. [5] 



avoid overcounting, the disallowed diagrams in this ex- 
pansion of E are the ones that can be cut into two by 
cutting two double lines in them. 

In order to have a better idea as to how this resum- 
mation happens, let us find all the single-link diagrams 
that contribute to the second-order self-consistent dia- 
gram of E in Fig [5] Note that this diagram comes with 
sign (—1). The double line between the cumulant of the 
two scatterings is built of single lines and cumulants of 
the scatterings between them. Expanding the double line 
in terms of single lines, diagrams with m scatterings have 
sign (— l) m in G. Therefore their contribution to the in- 
vestigated diagram has sign (— 1)(— l) m . 

In the single-line expansion of E lets consider a subset 
of diagrams within the (m + 2) nd order scatterings. Let 
this specific subset consist of diagrams which have the 
second-order cumulant of the first and the last scattering 
(i.e. the crosses representing them are connected to each 
other by a dashed line but not to any other crosses) and 
all the possible combination of cumulants of the inter- 
mediate m scattering potentials. Note that these are all 
parts of the self-energy since they cannot be cut to two 
by cutting one of the single lines in them. The sign of 
these diagrams is (— l)( m + 1 ) which agrees with the sign of 
the m-th order term of G above. Therefore, we found all 
the single line diagrams that contribute to the mth order 
term of G in the self-consistent diagram investigated. 

The above procedure of replacing single lines with dou- 
ble lines is called self- consistent perturbation expansion 
in the literature of field-theory. Assuming that the ex- 
pansion works, to obtain the n th order correction to the 
self-energy one has to calculate the Green's function us- 
ing the (rt — l) st order approximation of E then substitute 
this GF into the nth order diagram of E. In other words, 
E has to be calculated using successive approximation. 

As done in the single-line case, let us compare the 
first- and second-order diagrams of the self-energy, de- 
noted by E 1 and E 2 , to check the validity of this new 
self-consistent expansion. As long as the higher-order di- 
agrams are higher order in the expansion parameter (p 
or q) the perturbation expansion works. 

Without loss of generality, let us investigate the long- 
distance behavior of two diagrams: the first-order dia- 
gram, the annealed diffusion operator, 



£«(r) 



qp 

Nr a 



(63) 
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and the second-order diagram, from Eq. 



£( 2 >(r) 



2(G(0) - G(r))q 2 p 
Afr a 



(64) 



Note that, in £ 2 (r), (G(0)—G(r)) is approximated to first 
order, i.e. it is calculated with only the first diagram in 
the self-energy. The details of its calculations are given 
in Appendix [Cl 

In one dimension, there are four different cases: 
1) a < d; (G(0) - G(r)) oc (pg)-* and 



S 2 (r)(x-W W 



(65) 



X 2 (r) <C S 1 (r) if q <C p. The perturbation expansion can 
be applied in the weak interaction limit, q — > 0, keeping 
the number of links fixed, but breaks down in a network 
with a few but strong long-range links, as p — > 0. 
2) d< a < d + l: (G(0) - G(r)) oc (pg)^ and 



>pq 



(66) 



S 2 (r) <C S x (r) if q 2 ~ a << p. A similar condition as in 
case 1) though, here, the strength of the interaction must 
have a smaller value for the expansion to be applicable. 
3) d+l<a<d + 2: (G(0) - G(r)) oc j^r a ~ 2 and 



£ 2 (r) oc f 2 . 



(67) 



The perturbation expansion breaks down for all param- 
eter values. 

4) d + 2 < a: The perturbation expansion breaks down. 
Though, unlike in the cases above, in this regime the ef- 
fect of the random links is negligible at large distances 
to that of the regular diffusion on the chain in both the 
annealed and the quenched system. 
In two dimensions there are two cases: 

1) a < d + 2: (G(O)-G(r)) oc -\n(pqa 2 ). Therefore, the 
condition for the validity of the perturbation expansion is 
q~~ 1 e~ 1 / q <C p. The perturbation expansion works when 
q — ► 0. In the limit of a small number of strong links, 
the situation is better than in the one-dimensional cases 
since the expansion is valid for a large range of p values, 
but still breaks down in the p — > limit. 

2) d + 2 < a: Similar to case 4) in one dimension, the 
perturbation expansion breaks down, though the diffu- 
sion through the random links are irrelevant at large dis- 
tances compared to the regular lattice diffusion. 

In three and higher dimensions (G(0) — G(r)) oc a d ~ 2 
and 



£ 2 (r) oc 



2 r\— 2 

pq z a z 



(68) 



The expansion can be applied in the weak interaction 
limit. In the case of vanishingly small number of links, 
unlike in the lower-dimensional case, the higher-order 
terms do not dominate the behavior of the self-energy. 



In fact, all the higher-order terms are the same or higher 
order in p as £ . The appearance of higher-order dia- 
grams with same p-dependence is due to the fact that the 
cumulants of A rnd are all the same order in p. Though, 
diagrams with multiple cumulants (see Fig. [9]) are of or- 
der p 2 . As it will be shown later, the approximation of 
the first order (annealed) diagram is valid in this case, 
the perturbation expansion provides the right scaling of 
the self-energy. 

In lower dimensions, d < 2, the failure of the self- 
consistent perturbation expansion in terms of the expan- 
sion parameter p is due to the fact that diagrams with 
multiple scatterings off a single link (i.e. diagrams with 
only one cumulant in them) contribute with the same 
weight, p, to the self-energy. These scattering processes 
also contain the GF propagating between two scattering 
events. The GF is determined from lower order approxi- 
mations of the self-energy, and therefore is a function of 
p itself. In the p — ► limit, the GF should converge to 
the original pure GF, G°, that is divergent as a function 
of the system size in these dimensions. Therefore, G is 
divergent as well when p — > 0. Following this argument, 
one can understand that the higher-order single-link (or 
one-cumulant) diagrams dominate the small-p limit be- 
havior of the self-energy since the first-order, single-cross, 
diagram does not depend on the GF, while all the other 
single-link diagrams do. 

As it was done in Subsection IIII Al the infinite set of 
single-link diagrams [85j can be summed. Using Eq. (|59p . 

oo 

7i—l a<b 

(2q(G aa - Ggb^qp 
X M]a-b\ a 

= - 5ib)(8gj - Sbj) 



a<b 



Af\a-b\«l + 2q(G aa -G ab )' 



(69) 



Since the operator is translationally invariant and the 
continuum limit is meaningful one can write 



E si (r) 



-P 



Nr a 1 + 2g(G(0) - G(r)) 



(70) 



Since the self-energy is a conservative (diffusion) opera- 
tor, 



£ s( (0) = -^E si (r). 



(71) 



This property can also be derived from Eq. (f6"9"]) . The 
above pair of equations is referred to as the Self- 
Consistent Formula (SCF). 

In order to solve the self-consistent formula, Eq. (|70p . 
the following steps are to be taken: 

First, an ansatz has to be made about the large-distance 
behavior of the self-energy. In most cases, we can make 
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a very good guess by assuming that the long-distance 
spatial behavior of the self-energy will be the same as 
that of the self-energy of the annealed system, 

£ s '(r) = -1 oc £ an (r) = [A™ d ](r) if |r| > a, (72) 

and the strength of the interaction, s(p, q), will be some 
function of p and q determined by the self-consistent for- 
mula. 

Second, (G(0) — G(r)) has to be calculated using the 
above ansatz, G(r) = (A + S si ) _1 (r). As done before, 
in order to calculate G(r), the Fourier transforms of the 
operators are used because of their diagonal form in k- 
space: 



G(r) 



1/a 



d d k 



l/L 



S(fc) (2ir) d 



and 



(G(O)-G(r)) 



i/». 



1 -e 



kr 



d d k 



1/L fc 2 + E(fe) (2tt)<* 



(73) 



(74) 



The calculation of the scaling of the latter quantity can 
be found in Appendix [Q for various forms of E(fc). 
Third, the asymptotic form of (G(0) — G(r)) is substi- 
tuted into the self-consistent formula 86]. If our ansatz 
about £(r) was correct in the first step, the self-consistent 
formula has to be satisfiable by choosing s appropriately; 
if it was not correct, a new ansatz has to be made. 
Fourth, having S s ', G(0) is calculated to obtain physi- 
cally measurable quantities of the system. The details of 
the calculations are as follows: 



E. The one-dimensional phases 

The leading-order results for one dimension are 
sketched in Fig. [21 Similar phases appear as were pre- 
dicted by the annealed argument. But, in the tran- 
sient/smooth phases, we have a different scaling prop- 
erty - a faster divergence as p vanishes. Furthermore, the 
recurrent/rough phases have different phase-boundaries. 
Recurrent/rough phase I, where the effect of the long- 
range links is negligible, spans a wider interval of the 
a-axis, and recurrent/rough phase II is collapsed to one 
point on this axis. In this phase, the sublinear behav- 
ior of G(0) is no longer determined by the distance- 
distribution of the random links but rather by their den- 
sity, p, through Eq. (fi?2")) . The form of this equation is 
approximate, but the fact that the exponent fi depends 
continuously on p is likely to be exact, in light of the 
scaling arguments in Section [IVI 

For < a < 1, let us go through the process step 
by step. The ansatz, Eq. ([72")l . is S(r) = and its 
Fourier transform is S(fc) « s. In the second step, from 
Appendix E2 



(G(0) - G(r)) a 



1 



(75) 



for large r-s. Substituting this asymptotic form into the 
SCF, Eq. CZ0]), 



P 



Afr c 



Therefore 



1+3 



s oc p 



s-)-0 P\fs 



(76) 



(77) 



In the last step we assumed that s — * as p — > 0. This 
assumption, which is indeed satisfied, can be checked for 
consistency after obtaining the functional form s(p). In 
the final step, from Appendix [Bl 



G(0) 



i/« 



l/L 



dk 



- p A 



oc p 



(78) 



For 3 < a, the Fourier transform of the ansatz is oc 
k 2 . Therefore, 



(G(0) - G(r)) oc r. 
Substituting this in the SCF, 



P 



1 



Afr a Nr a 



(79) 



(80) 



The spatial behavior of the two sides of the equation are 
different, our ansatz does not work. This problem can be 



fixed easily by setting £(r) 



The Fourier trans- 



form of any self energy that decays faster than r ( d + 2 ) is 
£(fc) oc k 2 , therefore (G(O) - G(r)) oc r is still valid for 
the new ansatz. Using this new form of E(r), the self- 
consistent equation is satisfiable if s oc p. The scaling 
behavior of the GF is G(0) oc L. 

For 1 < a < 2, the Fourier transform of the ansatz is 
£(fc) oc sk - 1 . Therefore, 



(G(0) - G(r)) oc s —° 



(81) 



Substituting this asymptotic form into the self-consistent 

3-a 

formula, s oc pi-<* and 



G(0) oc 



(82) 



For 2 < a < 3, the Fourier transform of the ansatz is 
oc s/c Q_1 , which is the same as in the previous case 
but G(r) has a different long-distance property: 



(G(O)-G(r)) oc -r a - 2 . 

s 

Using this result in the self-consistent formula, 

s , P s 
Afr a Mr a r a ~ 2 ' 



(83) 



(84) 



A better choice is S(r) oc ^£pr- Since a + 1 > 3, its 
Fourier transform is oc sk , and therefore (G(0) — 
G(r)) oc r. The SCF is satisfiable if s oc p, resulting in 



G(0) oc L, 



(85) 
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For a = 2, the Fourier transform of the ansatz is 
k. Therefore, 



(G(O)-G(r)) cx - ln(rs) 



Substituting into the SCF, 



s ^ P s 
J\fr Mr ln(rs) ' 



(86) 



(87) 



the two sides have a different r dependence. One iteration 
step generated a faster decaying function. This tells us 
that the fixed point of the SCF has to be a function 
with a faster decay than that of out initial ansatz at long 
distances. Let us modify the ansatz of E, to £(r) oc -f-^, 
where 2 < 7 < 3 and n is a constant factor. For self- 
energies with such an exponent, (G(0) — G(r)) has already 
been calculated in the 2 < a < 3 case above. As it will be 
shown below, in this phase the (7 dependent) constants in 
E(r) will have a crucial role in determining the exponent 
7- 

For this reason one has to calculate the scaling prop- 
erties of the ansatz of the self-energy with more scrutiny: 
from Appendix [Althe result of the Fourier transformation 



E(fc) = — k~<- 
ni 



where n\ = n (2sin(§7)r(l-7)) . Also, from Ap- 
pendix [Cj using the above form of the self-energy, 



(G(0)-G(r))«^ 
s 



(89) 



for large distances, where 



n 2 



n cos (§7) T(2 - 7) 
2tt sin(§7)r(l -7) 



71 / 7T > 

2^ COt (r) (1 



-7) 



Substituting this into the SCF, 



s 
nr 1 



P 



s p s 1 

Mr 2 n 2 r~<~ 2 Mn 2 r^ 



(90) 



(91) 



The two sides of the equation have the same functional 
form. Though s is undetermined by the single link ap- 
proximation, we will see later that some information is 
still known. The only remaining tunable parameter is 7 
(in n 2 ), to equate the two sides of the formula. Noting 
that M — 2J2ili = 'Xj the equation determining the 
value of 7 as a function of p, 



l=p 



3 tan(f 7 ) 
2?r 7-I 



, where 2 < 7 < 3 



(92) 



has to be solved for different p-s numerically. Therefore, 
G(0) cx -L^p)- 2 . (93) 



Note that even though the single-link approximation 
does not give information about s(p, q), one can still ar- 
gue about its q dependence: The dimension of s has to 
be that of (length) 1 " 3 since k 2 and sfc 7-1 should have 
the same dimensions in G(k). Since p is dimensionless 
for d = 1 and a — 2, q with dimension (length)" 1 is 
the only parameter of the problem that can restore the 
dimensionality of E(fc). Therefore one may expect that 
s(p,q) oc q 3 ~ 7 .f(p) where f(p) is undetermined. 



F. The two-dimensional phases 

The results for the quenched two-dimensional system 
are summarized in Fig. [3J The overall difference between 
the annealed and the quenched cases is that the quenched 
self-energy acquires logarithmic corrections to that of the 
annealed system. As a result the quenched phases have 
stronger (or equal) logarithmic divergences than their an- 
nealed counterpart. Three different phases are identified: 
a logarithmically smooth/transient phase when a < 4, a 
sub-logarithmically rough/recurrent phase when a = 4, 
and a logarithmically rough/recurrent phase when a > 4. 
To map the phase space, the single-link self-consistent 
formula [Eqn. ({TO")) ] is used to obtain the leading order 
behavior of the GF. 

The steps of solving the self-consistent formula are the 
same as in the one-dimensional case. The details are as 
follows: 

For a < 2 the Fourier transform of the ansatz is E(fc) oc s. 
Therefore, 



(G(0) - G(r)) oc ln(sa 2 



(94) 



and the result of the self-consistent equation as s — > is 
P 



so 



and 



In s 



bxp 



G(0) oc In 



(lnp) 



(95) 
(96) 

(97) 



For 2 < a < 4 the Fourier transform of the ansatz is 



T,(k) (x sk a ~ 2 . Therefore, 



(G(0) - G(r)) oc \n(sa 



4— a\ 



(98) 



and the result of the self-consistent equation is the same 

as in the previous case: s w + and G(0) oc In (^j^j ■ 

For a > 4 the Fourier transform of the ansatz is T,(k) oc 
sk 2 . Therefore, 



(G(O)-G(r)) 



1 - Jo (At) kdk 



1/L (l + s)k 2 (2n) 2 
f 1 / - dk 

\ -=Hr/a). (99) 

Jl/r k 



1G 



Substituting this to the self-consistent formula, 

s i P 1 
Afr a Afr a \n(r/a) ' 



(100) 



Since the spatial behavior of the two sides of the equa- 
tion are different, our ansatz does not work. As in the 
one-dimensional problem, this problem can be fixed by 
setting S(r) = jf^r^rj^- The Fourier transform of any 
self-energy that decays faster than y-( d + 2 ) is £(fc) oc k 2 , 
and therefore (G(0) — G(r)) cx ln(r/a) is still valid for 
the new ansatz. Using this new form of £(r), the self- 
consistent equation is satisfied if s ~ p. The scaling be- 
havior of the GF is 



G(0) cx InL. 



(101) 



For a — 4, let us look into the details. In real 
space, E(r) = jf^, and its Fourier transform is S(fc) = 
sk 2 ln(fca) (see the calculations of Appendix |A")) . There- 
fore, from Appendix ICl 



(G(0) - G(r)) cx - ln(s ln(r/a)). 



(102) 



Substituting the above result into the self-consistent 
equation, 



V 



r 4 r Afr 4 ln(sln(r/a))' 



(103) 



The two sides of the equation do not match because 
the iteration process generated a double-logarithmically 
stronger decay of the self-energy. If one tried another 
ansatz, S(r) cx e <C 1, the iteration step would 

weaken the decay making it cx r 4 i n ( r / a ) ■ Therefore the 
fixed point of the self-consistent equation for a = 4 
should be a function with a decay between ^ and -^hr;. 
After some guess work, a better ansatz is 



E(r) 



n r 4 hr" (r/a) ' 



(104) 



where > fi > 1 and n is a constant. Its Fourier space 
behavior from Appendix IA1 for small k-s is 



£(fc) 



SCi(/x) 



fc* In 1- " (ifea) +0 



k 



(in" 



(105) 



where ci(/i) is a constant, the result of the Fourier trans- 
formation. From Appendix [Cl 

G(0) - G(r) oc ^4^T ln M (r/a) + O (in"- 1 (r/o)) , 
s ci(/x) 

(106) 

where C2(/i) is a constant, the result of the inverse- 
Fourier transformation. Substituting this result to the 
self-consistent formula when r — > oo, 



P 



n c 2 (m) 

X r-r X 



1 



nr 4 ln"(r/a) Afr 4 s a(ji) ln"(r/a) 



(107) 



yielding 



1= P_ x cijp) 
Af c 2 (m) 



(108) 



where A/" = ^ r ^ s * ne normalization constant of the 
probability distribution of the random links, and Ci(/x) 
and C2(p) are constants defined by the Fourier transfor- 
mation of specific functions. As in the one-dimensional 
case, fi became a continuously varying function of p, de- 
fined by Eq. (|108[) . Though, in this calculation, c\ and 
C2 are not determined, in future research one may try to 
perform the Fourier transformations more accurately or 
determine their value numerically by performing discrete- 
Fourier transformation. 
Finally, one finds 



G(0) 



l/a 



kdk 



\ /L k 2 + sk 2 \n 1 - f *(ka) 
oc iln" (p) (£), < n{p) < 1, 



(109) 
(110) 



a "sub-logarithmic" divergence with the system-size. 

From numerical and experimental point of view, one 
must note that the above logarithmic differences are dif- 
ficult to measure. 



G. Higher dimensions, d > 3 

In the naive perturbation expansion it was seen that 
higher-order single-link diagrams give the same order 
contribution in p as the first order one. From the SCF, 
one can verify this result since the GF is finite and deter- 
mined by its behavior at the lattice scale (see Appendix 



0, 



(G(0) - G(r)) » g cx a 



d-2 



(111) 



where g is some constant and the exact value of it is de- 
termined by the microscopic details of the lattice. Sub- 
stituting this into the SCF, 



E s '(r) 



P 



Afr a 1 + 2qg 



(112) 



and G(0) is independent of p and q to leading order, 
G(0)«5 + — 



H. Higher-order corrections 

In the previous sections the self-energy was approx- 
imated with the infinite sum of single-link diagrams. 
Here, the contribution of higher-order diagrams is inves- 
tigated in the case of the plain SW network (a = 0) and 
the other cases, a ^ 0, are left to future research. The 
first such diagram is that of Fig. [HJa), the lowest-order 
two-cumulant diagram. It is called a two-cumulant dia- 
gram since it has two sets of crosses: in each set, crosses 
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x x x x x 



x x x x x 



FIG. 9: Two cumulant diagrams. The notation is the same 
as in Fig. [5] Moreover, square brackets, [ ] were dropped from 
around double lines for the sake of clarity. 




FIG. 10: Sketch of the leading order processes in the two- 
cumulant diagram of Fig. |9ja). Links are represented by fat 
lines, GF-s are represented by wavy lines, and the underly- 
ing d-dimensional lattice is represented by a flat sheet. The 
two indices of the corresponding matrix are represented by 
straight lines with an arrow, a) the leading order process 
at short distances; (b) the leading order process at long dis- 
tances. 



are connected to each other by dashed line. According 
to the definition of the diagrams, these sets correspond 
to separate cumulants of the perturbation potential. 

In terms of matrices the corresponding scattering po- 
tential is 

E 2c = q 4 V" y(a.b) G y(e,f) G y(a,b) G y(eJ) 



a,b,e,f 



(113) 



sake of simplicity and [(x^ e, ^) n ] = p/L d for plain SW 
networks. Note that each operator in the sum is con- 
servative [see Eq. (f57j) ] hence the whole diagram is so 
too. Let us investigate a single term of this sum. Using 
Eq. (|48]), the factorized form of V^'^ and V^ eJ \ one 
can see that such a term consist of three scatterings of 
the GF back and forth between two links in the system: 
a link between a and b and an other one between e and 
/. All the possible combinations of scatterings between 
the links are present. It can be shown and we will see it 
in an example that the most important terms in the sum 
are such where all the three GF-s in the diagram scatter 
between the same legs of the links, let us say between a 
and e. At small distances, the most important process 
of them is when the process starts with scattering off a 
then the GF scatters between a and e three times and 
finally it "scatters out" at e (see Fig.fTOla)). The matrix 
corresponding to this process is 



a,e b,f 



(114) 



At large distances the main contribution comes from a 
process depicted in Fig. HW b) [S7| ] . It starts off exactly 
the same as in the previous case but in the final step the 
particle jumps over the long-range link between e and / 
and scatters out at /. The corresponding matrix is 



^2c,B 



-^EE<u G -) 3 ^ 



a,f b.e 



aj b e 

-? 4 f^E( G -) 3 - ( 115 ) 



Note that these two processes are not conservative, yet 
the sum of them is, which is required for the self-energy. 
In summary, in the continuum limit, 



1 fB (r) = qV ( G{rf - j- d I G(r') 3 dV 



Note that for long distances 

S 2C (r) « -<? 4 ^ / Gir'y'W'r' . 



(116) 



(117) 



where the indices of the matrices were dropped for the 



the diagram is distance independent and uniform, and 
the contribution of it is a mean-field interaction just like 
the contribution of the leading order terms but with a 
different strength. 
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b) 



FIG. 11: Sketch of a process the contribution of which van- 
ishes in the thermodynamic limit compared to the leading 
order one, depicted in Fig. [TU] 



As an example, let us investigate another term in the 
diagram, where not all the GF scatter between the same 
legs of the two links, and see that the contribution of 
such a diagram is vanishing in the large system size limit. 
This process starts with scattering off a. The GF then 
propagates back and forth between a and e twice before 
it jumps over the long-range link between a and &, where 
it then propagates from b to /, and finally jumps from 
/ to e (see Fig. [TTJa)). The corresponding matrix, up to 
combinatorial factors, is 



b f 



bf 



(118) 



/ 



and the long-distance counter-term of it is [see Fig. fTTT b)] 



f b 



(119) 



In the continuum limit 



(r) = q 4 ^[ (CM 2 / G(r')dV 



G(r') 2 dV / G(r")dV 



(120) 

Compared to the leading order term, H 2c,A+B (r) , this 
matrix vanishes in the thermodynamic limit. 



FIG. 12: Resummed two-cumulant diagrams. Stars represent 
the resummed scattering potentials introduced in Sect. IIIlTTl 
Otherwise, the notation is the same as in Fig. [5] Square 
brackets, [ ] were dropped from around double lines for the 
sake of clarity. 



Qualitatively, the above result can be understood from 
the following argument: For the process depicted in 
Fig. [TUJ only two endpoints, a and e, of the two links 
have to be within a distance proportional to the corre- 
lation length of G of each other so the GF can propa- 
gate with finite weight between them. For the process of 
Fig. [TTJ both endpoints of both links have to be within 
a distance of the correlation length from each other, a 
close to e and b close to /. Since the two links are ran- 
domly spread on the lattice in a uniform fashion, having 
both endpoints of both links in the vicinity of each other 
is much less likely than having only one endpoint of one 
link close to one endpoint of the other link. 

In low dimensions, some higher-order two-link dia- 
grams have the same problem as the higher-order single- 
link diagrams. For single-link diagrams, it was demon- 
strated that the second-order diagram was more domi- 
nant than the first-order one in the p — ► limit. There 
are two-link (or two-cumulant) diagrams that involve 
GF-s propagating between the endpoints of the same link, 
just like in the case of single-link diagrams. For exam- 
ple, the last double line of Fig [HKb) represents such an 
event (it can be checked by expressing the diagram in 
terms of matrices). Such a GF will have a contribution 
to the diagram which is divergent as p — ► 0. Since this 
divergent contribution is the only difference between this 
diagram and the diagram in Fig. [9](a), the higher-order 
diagrams will dominate the small p behavior. As it was 
done before, this problem can be avoided by resumming 
over all the single-link scatterings of such diagrams. As a 
result, when there is a scattering off a link, the scattering 
is substituted by all the single link scatterings. In terms 
of matrices the substitution 



g I/ (ab) -> qV {ab) - qV {ab) GqV {ab) 
+ qV {ab) GqV {ab) GqV {ab) - 



(121) 



is made for all x-es in the diagrams. It is known from the 
single-link calculations that the above infinite sum yields 



y(ab) 



1 + 2q(G aa - G ab ) 



(122) 



One can say that the result of the resummation is that 
the strength of the interaction gets renormalized from q 



to 



l+2q(G aa -G ab )- 

In a typical situation the separation of a and b diverges 
with the system size. Therefore, in such situations, G a b is 
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© = [*]- 



I* + *z 



FIG. 13: The self-energy in terms of the resummed potentials. 
The notation is the same as in the previous figure. 



vanishingly small because the GF has a finite correlation 
length. As a result, the resummed scattering potential 
can be well approximated by 



y(ab) 



1 + 2qG{0) 



(123) 



in the continuum limit. In diagrammatic notation, 
the above resummation is represented by replacing the 
crosses with stars as shown in Fig. Q21 The expansion 
of the self-energy in terms of the resummed potentials is 
shown in Fig. [T3j 

Using these renormalized single-link scattering poten- 
tial, let us calculate the leading order behavior of the 
diagram in Fig. H27 a). The calculation is the same as 
was done for the diagram of Fig. [U[a) except that q is 
replaced by 1+2 %c(o) in Q 116 ! 1 • 

In one dimension, G(0) oc p -1 . In the limit p 0, 



<1 



1 + 2qG{Q) 



P- 



(124) 



Therefore, the second-order term in the resummed self- 
energy scales as 



£ 2c (r) oc - 



A 9 
P P 



G{r'f&r' 



(125) 



at large distances. The integral in the expression can be 
approximated using the results of Appendix [Dj 



s -l/2 tf r< £ 



G(r) cx 



(126) 







if r>£ 



where s oc p 2 , obtained from the single- link approxima- 
tion [see Eq. ([77| ], and £ = s -1 / 2 is the correlation length 
as before. Therefore 



and 



G{r'fAr' oc p~ 



E 2c (r) oc 



(127) 



(128) 



The Fourier transform of this operator for small 
wavenumbers £ 2c (fc) oc p 2 has the same p dependence 
as the leading order one, Eq. (|TT|) . Higher-order dia- 
grams, like the one in Fig. [T2Tb). yield contributions of 
the same order, p 2 . As a result, it can be concluded 
that the self-energy scales as p 2 for small p-s, though nei- 
ther the coefficient of it nor higher-order corrections are 
known. 



In two dimensions, G(0) oc lnp from Eq. (|9"7) , therefore 
q is renormalized to 1/ lnp. The second-order correction 
to the self-energy at large distances 



£ 2c (r) oc 



G( 



~') 3 dV. 



(129) 



(lnp) 4 L 2 

From Appendix [Dl the integral can be approximated us- 



ing 



const. + ln(r/£) if r -C £ 
Gin x { (130) 
if r > £ 



where S, — s 1 / 2 oc y/\np/p from Eq. (j9"B"|) . Therefore, 



G(r')Vdr' oc 



lnp 



and 



S 2c (r) oc 



as p 



Qnp) 3 L 2 ' 



(131) 



(132) 



0, S 2c (fc) oc p/(lnp) 3 < 



In Fourier space, 
S s '(fc) oc p/lnp. 

In conclusion, in two dimensions, the single-link per- 
turbation expansion of Fig. [T3] works. The mean-field 
model acquires logarithmic corrections from sample-to- 
sample fluctuations. The fact that the perturbation ex- 
pansion works in two dimensions raises the possibility of 
an e-expansion for d = 2— e which was studied in [23| . For 
small e-s, it was shown that the higher-order diagrams of 
the perturbation expansion produce higher-order terms 
in e. Though e = 1 for d = 1, the e expansion still defines 
an ordering of the diagrams making it possible to obtain 
more accurate analytic predictions for the behavior of the 
system in the p — > limit. 



IV. SCALING ARGUMENTS 

After the more careful perturbative treatment of the 
previous section, let us revisit the effect of long range 
links and estimate the behavior of the quenched system 
using scaling and renormalization arguments. 



A. Estimating the validity of the annealed 
approximation 

In many statistical systems introducing random long- 
range links to the underlying lattice results in mean- 
field (MF) or, in our terminology, annealed behavior [381 ] . 
In the following, we will study the self-consistency of 
the annealed approximation by estimating the samplc- 
to-sample fluctuations of the quenched systems around 
it[z| 

First, we consider the plain SW network (a = 0). For 
this estimation, we utilize the EW Hamiltonian [Eq. 
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expressed in terms of the random component of the ad- 
jacency matrix, A Ttr i , and the difference of the fields at 
site r and r': 

H[h] = H a [h]+q^2A r . r ,(h T - h r , f 
= H [h] + q^2[A r . T ,}(h r - h r ,) 2 

v.r' 

+ qJ2(A r ^ - [A r ^})(h r - h r ,) 2 

r.r' 

= Ho [h] + [H rnd ] [h] + H™tt. [h] , (133) 

where Hq is the Hamiltonian of the EW model on the 
regular, unperturbed, lattice and H rnd is the "energy" 
due to the random part of the network. 

According to MF theory, the last term in H, 



(134) 



is negligible compared to the second one because the 
sample-to-sample fluctuations of A r y are suppressed due 
to averaging over a large number of sites. For a link be- 
tween two sites, the relative size of these fluctuations can 
be quantified by 



y/[(A ry )*} - [A r yY< 
[A r y] 



(135) 



the ratio between the standard deviation and the mean 
of the random potential. 

For the whole system, one has to average over many 
sites to estimate the magnitude of the fluctuations in the 
energy, H[h\. Averaging over the entire system would 
suppress the fluctuations, though doing so would under- 
estimate their strength. This can be understood in the 
following argument: 

If the separation of r and r' is less than the correlation 
length, £, of the field h r , the values of the field at those 
two sites are close to each other, and (h T — h r /) 2 is small. 
Therefore, most of the contribution to the energy comes 
from sites where r' is outside of the correlation volume 
of the field at site r. 

On the other hand, the fluctuations in the energy will 
only amplify the effect of each other where the values of 
the field are correlated, i.e. over volumes oc £ d . There- 
fore, the condition for the MF or annealed approximation 
to be valid is 



/ d d r J dV([(V) 2 ]-Krf) 

|r|<£ £<|r-r'| 

/ d d r ] d d r'[A T y] 

|r|<£ £<|r-r'| 



< 1. (136) 



Here, we used the fact that the standard deviation of 
a sum of random variables is the square root of the 
sum of their variances. From Subsection IIII C( L4 r . r '] = 
([(^4r,r') 2 ] — L4 r ,r'] 2 ) = P / L d and the correlation length 



is £ cx (pq) . The above condition can be expressed 
as 



P 



-X+d/2 d/2 ^ 1 



(137) 



The assumptions of the MF approximation are always 
valid when q — ► 0. In one dimension, if p <C q, the 
sample-to-sample fluctuations destroy the MF behavior 
otherwise MF approximation can still be applied. In two 
dimensions, if p — *■ the MF predictions acquire logarith- 
mic corrections from the sample-to-sample fluctuations 
(see Section fill Fj) . In three dimensions and above, the 
MF criterion is always satisfied. 

Second, when a < d, [H rnd ] generates a mass (or finite 
correlation length) just like in the case of a plain SW 
network and £ cx (pq) as well. The same argument 
can be repeated as in the previous case. The condition 
of the MF behavior, Eq. fBS]) . is 



p -l+d/2 q d/2 <<; L 



(138) 



Third, when d < a < d + 2, strictly speaking, there 
is no correlation length in the system, though the typi- 
cal size of correlations is proportional to l x , defined by 
the lengthscale above which the long-range random links 
dominate the behavior of the system. In Fourier space, 



(139) 



where k x = l x . Thus, for the lengthscale l x one finds 

l x = (pq)~ 2 + d -° . Repeating the same argument as 
above, replacing the role of £ by l x , and using [j4 rr /] = 
([(^4 r ,r') 2 ] — [A-.r'] 2 ) = p/(A/V Q ), one arrives at the con- 
dition 



P 



-2 q 2d- 



< 1 



(140) 



for the MF approximation to be valid. Observe that in 
one dimension, if 2 < a < 3, the MF approximation 
breaks down for all p-s and q-s. In all the other cases the 
condition is similar to those of the plain SW network, 
but due to the distance-dependent cutoff of the links, a 
higher density of long-range links is required for the MF 
approximation to be valid (i.e. p has to be larger for 
observing MF behavior) . 

Fourth, when d + 2 < a, the typical size of the correla- 
tions is that of the system size, so one cannot apply the 
above argument. Though, one can introduce a "virtual" 
correlation length to the system, £ v , check the condition 
of the MF approximation, Eq. (|136[) . as a function of 
and in the final step take the limit when the correlation 
length goes to infinity (or, at least, to the system-size). 
The result is 



P 



-lt(« 

sv 



-2d) 



< 1 



(141) 



In conclusion: For d + 2 < a < 2d, the condition for the 

MF approximation is satisfied since tii a - 2d) - as Cv - 
oc. For 2d < a, the condition of for the MF condition is 

not satisfied since £y" 2d ^ — > oo as £v - * oo. 
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FIG. 14: The steps of the renormalization procedure: 1) 
Take the original system. 2) Divide it into boxes of size I. 3) 
Average over the behavior within each i-sized box. 4) Rescale 
the system so that I — > a. 



B. Renormalization 

As it was shown above, MF theory is insufficient to de- 
scribe diffusive processes on quenched SW networks for 
certain parameter regimes. Though, renormalization ap- 
proach can capture the universal behavior of these sys- 
tems: Assuming that the large scale properties of the 
system are not influenced by the microscopic details, one 
can average over the small lengthscale behavior according 
to the steps of Fig. [l~4l 

In the following, we will investigate how the diffusion 
operator changes doing the steps of renormalization. Dif- 
fusion on PL-SW networks is described only by two pa- 
rameters: p and q (besides d, the dimensionality of the 
system). First, the change of the density of random links 
under a change in scale is considered. Second, we will 
look at the renormalization of q for only a single link in 
the system. Third, we will combine the results at hand 
to calculate the scaling of the GF on a PL-SW network; 
since such a network includes many links, we will have 
to include interaction between links when the density of 
links becomes of order unity at the given scale. 



where the normalization constant is 

i f L d-ij { L d - a ifa<d 
[ a d - a if d < a 

Therefore, 

( l d iia<d 
p oc p x < (144) 
( l 2d - a if d < a. 

In the fourth step, the (al) d sized cell is rescaled to size 
a d and the probability of one such node to have a random 
link will be that of the original cell, (al) d . In summary, 
the flow equation of the parameter of the probability dis- 
tribution of the random links is 

{l d \la<d 
(145) 
pd-a if d < a. 



D. Renormalization of q for a single sink and a 
single link 

In renormalization group theory, the procedure to ob- 
tain the universal long-wavclcngth properties of the sys- 
tem is averaging over short wavelength modes. Due to 
interactions which mix long and short wavelength modes, 
the Hamiltonian that describes the statistical properties 
of the system will change: these mixing interactions will 
generate modified or new interactions between long wave- 
length modes after averaging. 

Let us follow through this procedure in the case of 
when there is only a single sink in the system. The 
Hamilton operator of the system is H = H°+qV. The un- 
perturbed Hamiltonian is the diffusion operator, H® k , = 
Skk'k 2 in Fourier space. The perturbation potential is 
qVij = qdiodoj in real space and 



C. Renormalization of p 

As defined in the Introduction, p is the probability of 
a site having a random link emanating out of that site 
irrespective of where the link ends. The first two steps 
of the renormalization procedure are straightforward. In 
the third step, after averaging over a cell of size (al) d the 
probability of a node within the cell to have a random 

link pointing outside of that cell is p r \f r dr ■ For all 
the nodes within the cell, this probability is 



r d 1 dr 



oc I p X 



"AT 



(al)" 



TT 



if a < d 



if d< 



(142) 



L d 



(146) 



a "uniform" matrix in Fourier space. We will distinguish 
between long and short wavelength modes 



if k < l/l 



(147) 



if k > l/l. 



Making this distinction, let us introduce the following 
division of V in Fourier basis: 



V = 




(148) 
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where V\ and V2 act on the long and short wavelength 
modes and Va and Vb are mixing between them. The 
Laplace operator is divided likewise: 




(149) 



The usual method of doing the averaging is by consider- 
ing the partition function: 



X / £>$" e *"( H '2+9Va)#"+*'(gVB)*"+*"(gVk)*' < ( 15 o) 



After averaging over the short wavelength modes by per- 
forming the integral J 2?$" (...), a <£>' dependent expres- 
sion will remain which is the generated interaction due 
to the short wavelength "mixing" . The result is a new 
renormalized perturbation potential, q r enV ren , defined 
by 



Z = X>$ e * (H+ « v)<I> = / 'e* ' {Hl 



(151) 

We will see that for the single sink problem q re nV ren has 
the same form as V\ only the interaction strength gets 
renormalized q — > q ren . The above mentioned procedure 
of performing the short-wavelength integral could be car- 
ried out but it was merely introduced to recapitulate the 
basics of renormalization group (RG) transformation. In- 
stead, let us calculate q r enV ren using GF-s and the dia- 
grammatic technique introduced in Section IIII1 

Single lines represent G° divided into long-, G° , and 
short-wavelength, G° , propagators, represented by sin- 
gle lines with a slash and single lines with double slash. 
Single lines with a wave represent the renormalized GF, 
Qren _ _|_ q ren V ren )~ 1 . Crosses represents scatter- 
ings off qV and a cross in a circle represents q re nV ren 
(see Figs. H5] and HH). 

Qren -^[q^ [ s on iy defined for long wavelengths can 
be calculated using the original scattering potential, V. 
Scatterings of the long wavelength propagators can hap- 
pen in two ways: i) a long wavelength propagator can 
scatter off V to a long wavelength propagator [Fig.lTSTb)] , 
and ii) a long wavelength propagator can scatter off V to 
short wavelength propagators and then scatter "back" to 
a long wavelength propagator [Fig. [TBl c)]. i.e., the scat- 
terings of the long wavelength propagators are mediated 
by short wavelength propagators. 

In Fig. [TBFa), all the possible scatterings of the long 
wavelength propagators are represented. In Fig. 1Mb) 
the scatterings are reordered introducing a new symbol, 
a cross in a circle that represents all the possible G° me- 
diated scatterings of G° -s. Figure fTtJTc) is the reorder- 
ing of FigfTBFb) to the similar form of a Dyson equation 
where one can identify ® as q re nV ren . 

The sum in the formula of <g) can be solved. Notice that 
this is exactly the same problem as that of Section IIIII 



d) 



FIG. 15: Introduction of the short and long wavelength 
diagrams. a) Decomposition of the GF to short and 
long-wavelength modes, b) Simple scattering of the long- 
wavelength GF. c) Scatterings of the long- wavelength GF me- 
diated by short- wavelength GF-s. d) A typical scattering with 
multiple scatterings between the long and short wavelength 
modes. 



except that in that problem the intervening propagator 
between scatterings was G° while here G° . Therefore 



l + 9 G "(r = 0) 



Vl = q r enV T 



(152) 



where 



G°"(r = 0) = f 
Ji 



1/a d d k 




(153) 



The dimensionality of q is (length) d 2 [89j therefore after 
the fourth step of renormalization: 



q^q = q r enl (2 d) 
In one dimension, when I — > 00 , 

q 



1 + ql 



-I — * const. 



In two dimensions, 



5k irki = ^ 



For d > 2, 



q oc 



1 + qa 



2-d 



2-d 



<D(l 



2-d\ 



0. 



(154) 



(155) 



(156) 



(157) 



In conclusion, q is relevant in one dimension, 
marginally irrelevant in two dimensions, and irrelevant 
in higher dimensions. 

For a single-link, there are two cases to be differenti- 
ated: i) When the separation of the two endpoints (A) is 
much larger than the lengthscale over which the renor- 
malization is carried out (I <g; A). In this case, the vicin- 
ity of these endpoints can be considered as independent 
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FIG. 16: Obtaining the renormalized GF for the single sink 
and link problem, a) The perturbation expansion of the renor- 
malized GF, represented by a line with a wave, b) Introduc- 
tion of the renormalized scattering operator, <g>. c) The Dyson 
equation of the renormalized GF. 



rf-dimensional spaces for the short-wavelength GF to ex- 
ist in, connected by a link at the origin with strength q. 
As in the case of a single sink, the renormalized inter- 
action strength can be calculated using the diagrams of 
Fig. [TU As a result, 



1r 



9 



l + 2qG°"(r = 0) 



(158) 



and the RG equation for the interaction strength, q, is 
the same as for the single-sink problem: (Eqs. (|1551 11561 
II 5T|) ) . ii) When I 3> A, links and sinks get treated very 
differently: links can simply be ignored (q = 0) since 
the perturbation potential caused by a link has vanishing 
strength at these lengthscales while the effect of sinks can 
be still strong (see Eqs. (f!55l [T56l fT57)) ). 



E. The renormalization of q including interactions 
between links and the scaling of the Green's function 

We now combine the results for the RG flow of p in 
Eq. Q145|) with the results for q to describe the scaling of 
the Green's function on a small-world network. When p 
is small, we can use the single-link results for the scaling 
of q in Eqs. (|155ll56ll57p . When p becomes of order 
unity, however, these results will need to be modified. 

Consider a network with p and q both much smaller 
than 1. If a < 2d, then p is increasing under the flow 
and will eventually become of order unity. We will start 
by investigating this case. We restrict ourselves to the 
study of the case when d < a < 2d, which is the most 
complicated one, for the sake of brevity. In the omitted 
case, a < d, the derivation of the scaling properties can 
be obtained similarly using the results of the two previous 
sections. 

We begin by studying the case when d < a < 2d and 
d < 2 (see Fig. [T7]). As long as both p and q are much 
smaller than unity, q obeys the naive scaling q oc l 2 ~ d . 
This naive scaling for q will break down when either q 



becomes of order unity (due to multiple scattering off 
of a single link) or when p becomes of order unity (due 
to scattering off multiple links). So, we must determine 
whether it is p or q that becomes of order unity first. 
The link density, p becomes of order 1, from the naive 
scaling, at I ~ {l/p) 1 /^ 2d ~ a \ while q becomes of order 
unity, from the naive scaling, at I ~ {l/q) x /( 2 ~ d > . So, 
if (l/q) 1 ^ 2 -^ < (l/p) 1 / (2d - a '> then q becomes of order 
unity first. That is, if q S> p( 2 ~ d )/( 2d - Q ) then q becomes 
of order unity while p is still much smaller than unity. In 
this case, we can continue to use the single link results 
(|155[) for the RG flow of q and so q approaches a constant 
of order unity as the length scale continues to increase, 
while p still increases following Eq. (| 145[) . Thus, beyond 
this length scale, everything is set by the scaling of p and 
the original value of q becomes unimportant. Therefore, 
by dimensional analysis, using the scaling of G(0), we 
find that G(0) - p -(2-d)/(2d- a )^ Qn the Qther hand; if 

q <C p( 2 ~ d )/( 2rf-Q ), then p becomes of order unity while 
q is still much less than one. In this case, mean-field 
theory becomes accurate [9(j, and we find that G(0) ~ 

(p q )(d-2)/(2+d- a) ^ ThuSj 



q » p (2-«0/(2<*-a) ^ G(0 ) „ p -(2-d)/(2d-a) (15g) 



for 



As a consistency check, note that the two expressions 
G(0) agree when q - p (2-d)/(2d- a ) _ 

Next, consider the case that d < a < 2d and d > 2. For 
d = 2, we find that q is marginally irrelevant under the 
RG flow (|156[) . decreasing inversely with the logarithm 
of the length scale. For d < a < 2d, p will become of 
order unity at a scale I ~ {l/p) 1 /^ 2d ^ a \ The scale I is 
a power of p, hence the logarithm of the length scale is 
proportional to a logarithm of p. Thus, we can apply 
mean-field theory at this scale I with a j~ l/ln(l/p). 
Therefore, there are logarithmic corrections to mean-field 
theory for d = 2 and G(0) ~ ln(p/lnp). In the case 
when d < a < 2d and d > 2, q is irrelevant and we can 
apply naive scaling until the scale I ~ (l/p) 1 ^ 2 ^") at 
which point we can apply mean-field theory. Therefore, 
for d > 2, mean-field is always accurate in the limit of 
small p, q. 

Finally, there is the case a — 2d, when p is unchanging 
under the flow. For d > 2, q flows to zero, and then 
the effect of the links becomes negligible. However, for 
d < 2, q flows to a constant of order unity. Then we have 
a nontrivial fixed line of the RG flow with p arbitrary and 
q of order unity. This leads to a G(0) which varies as an 
anomalous power of the system size, where the power 
depends continuously on the value of p. 

Note that, when a > 2d, p always flows to zero there- 
fore the long-distance scaling of the GF becomes insen- 
sitive to the presence of the random links. 
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FIG. 17: The schematic representation of the RG flow when d < 2 and d < a < 2d for the quenched system. The parameter 
space is divided into two parts: i) where the behavior of the system is MF and ii) where the large-scale behavior is only 
determined by p. We also indicated the scaling of the GF in the two regimes. In part ii), consider two points in the phase 
space A and B with different q but the same p values: Renormalization transformed both systems to an identical one in the 
parameter space, where the GF is G*. Furthermore, it took the same number of renormalization steps, I, to get to G* from 
both points because I is only determined by p. As a consequence, considering that the GF can be well approximated by the 

free GF at small lengthscales (until p becomes of order unity): Ga(0) oc Gb (0) oc G*^ 2-d ' oc p 2d ~ a . In other words, G(0) have 
the same scaling behavior in the two points. 



V. APPLICATIONS AND NUMERICS 

A. Calculating physical properties from the 
Green's function 

In Section [IT] and IIII1 an impurity averaged perturba- 
tion expansion was set up to calculate the propagator of 
the diffusion operator for different configurations of the 
network by changing a, p, and q. Here we interpret the 
results for the GF in terms of the relevant observables of 
surface growth (Section II B[) and random- walk processes 
(Section HA)). 

For EW processes, from Eq. 1(25]). 

L d L d 

[<™ 2 >] = h ^ [Gu] = i?^ G(0) = G(0) - (160) 

i=l t=l 

Similarly in the annealed case [(w 2 )] ann = G ann (0). 
Therefore, all the results derived for G(0) can be trans- 
lated in terms of the average width of the EW process, 
represented by the naming of the different phases of Fig- 
ures [2] and [3] 

Though the perturbative calculations were done for 
G(0), they can be repeated for G(0,u>), defined in 
Eq. (|18p. The difference one has to consider is that, for 
the EW processes, the unperturbed GF is G° = (A ) -1 ; 
for random-walks, it is G°(oj) = (A + lu)~ 1 . This substi- 
tution also results in a slight change in the Fourier space 



calculations 

(161) 

as can be seen in Fig. where E(fc) is the self-energy 
generated by the random links as calculated in Sections ITT1 
and IIIII In order to obtain the scaling properties of 
G(r, u>) and, specifically, of G(0, ui), the approximations 
of Appendix IB1 are used. As a general rule, in the case of 
G(fc,tL>), the small-A; divergence is cut off by the u> term in 
the denominator. Using the approximation of Eq. (|19[) . 
the scaling of the expected number of returns can be re- 
lated to G(0,w= l/T) by 

[Fin cx ^]T[G r , r (a, = l/T)] + J 

r 

= ^£G(0,u; = l/T) + J 

r 

= G(Q, W = l/T) + -^. (162) 

Initially, we are only interested in the long-time behav- 
ior of the expected number of returns in the infinite 
system-size limit, therefore, [F(T)j oc G(0,u> = l/T) in 
the thermodynamic limit. As it was done for the width 
of the EW process, the interpretation of the scaling of 
G(0,u} = l/T) as the scaling of the expected number 
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of returns is represented in the naming of the different 
phases of Figures [2] and [3l Later in this section, we will 
also see how the finite size of the system changes the 
behavior of F(T). 



B. Numerical approach 

Numerical results were obtained in one-dimension 
(N = L). In order to check the analytic predictions of 
Section ILTT1 we numerically calculated different physical 
properties of the processes introduced in Section [T] For 
the quenched networks, numerical results were obtained 
by: (i) generating a large number of realizations (from 
100 to 1000) of PL-SW networks, then (ii) generating 
the diffusion operator, Eq. (UJ), on each realization, then 
(iii) using exact numerical diagonalization [79( to obtain 
the eigensystem of the diffusion operator on each real- 
ization, then (iv) using the spectral decomposition of the 
different physical quantities, Eqs. (|2U)) and (fTB)) , to ob- 
tain their value for each realization; then (v) averaging 
these quantities over a large number of realizations of the 
networks. 

For the annealed networks, we numerically integrated 
the corresponding timc-discretized stochastic EW pro- 
cess in a dynamically annealed network: at every time 
step the random links, with density p, are reassigned. 
The EW equation was numerically integrated using the 
most basic time-discretization scheme 

hi(t + At) = h i {t) + At[h i+1 (t) + h i - 1 (t)-2hi(t)} 
+ At qA i<r(i) (t) [h r{i) (t) - hi(t)] 
+ rn(t)V2M, (163) 

where unit lattice spacing for the underlying regular lat- 
tice is assumed, r)i(t)-s are independent and identically 
distributed random variables for all i and t with Gaus- 
sian distribution of zero mean and unit variance, and r(i) 
is the random neighbor of node i at time t. In this an- 
nealed construction, at every time step, each node has 
a random neighbor with probability p; Ai >r u\=X if node 
i has a random neighbor, A i r ^=0 otherwise. At every 
time step the random links are independently reassigned, 
a new PL-SW-network configuration is generated. For 
the time-discretization scheme to be convergent for the 
EW process on any fixed network, one can show that 



crete lattices), one has 



I At A n 



-II <1 



(164) 



i.e., At < 2/A max is required, where A max is the largest 
eigenvalue of the network Laplacian. While for the an- 
nealed network, where the network changes at every time 
step, we do not have a similar rigorous requirement, as 
a guiding estimate, one can consider the largest eigen- 
value of the average (annealed) Laplacian. For example, 
for a < 1 PL-SW networks, using our results for the 
annealed random Laplacian [Eq. (f3T))) ] (embedded in dis- 



A 



arm 
max 



max{2[l - cos(fc)] + [A rnd ]{k)} 



= max{2[l - cos(fc)] + Cqp} = 4 + Cqp, 



(165) 



where C is a constant which does not depend on p, q, 
and the system size N. The resulting time-discretization 
scheme then requires 



At < 



Cpq 



Similarly, for 1 < a < 3, one finds 



At < 



4 + Cpq4 a ~ 



(166) 



(167) 



The important consequence of the above order-of magni- 
tude estimates for the time discretization is that for an- 
nealed plain or PL-SW networks, At does not depend on 
the system size. Thus, if a sufficiently small At is tested 
and confirmed to lead to a converging scheme for a small 
system size, it will do so for all system sizes |9l| . In our 
numerical integrations, we used At = 0.10, which turned 
out to be sufficient to achieve numerical convergence for 
the values q = 1 and p < 1 considered here. 



C. Roughness scaling on plain SW networks 

In Section ITTT1 for plain SW networks (a = 0), it was 
shown that the MF approximation is valid in the high 
density of weak links limit (q — > 0) and [(u> 2 )] oc q^ 1 ^ 2 , 
but breaks down when the number of random links is 
small (p — ► 0). In the later case, the self-consistent per- 



turbation expansion yielded [(i 



oc p 



In Fig. [181 we 



compared our analytic predictions to numerical results 
and found a good agreement between the numerics and 
the analytics for an intermediate interval of q and p. For 
small values finite-size effects are observed, addressed in 
Section EU 



D. Roughness scaling in transient /smooth phase II 
and in recurrent/rough phases for annealed and 
quenched networks 

In Fig. [TH the different scaling behavior of the 
annealed and quenched system is shown in tran- 
sient/smooth phase II (a—lA). The quenched results 
were obtained by exact numerical diagonalization tech- 
niques, and the annealed results were produced by inte- 
grating the stochastic EW process in a dynamically an- 
nealed network, as described above. As can be seen in 
Fig. [ini the numerical results and the annealed analytic 
prediction are in good agreement. 

For the annealed case, Fig. [20] shows the system-size 
dependence of the width (or G(0)) in transient/smooth 
phase II (a=1.4) and in recurrent /rough phase II 
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FIG. 18: Disorder averaged width obtained by exact numer- 
ical diagonalization of the diffusion operator (a) changing the 
strength of the long-range links while keeping their density 
fixed (p = 1); and (b) changing the density of the long-range 
links while keeping the strength constant (q — 1) for system 
sizes indicated in the figure. The two slopes indicate the ana- 
lytic predictions in the two cases. Note that the MF argument 
would predict the same behavior for both cases. 

(a=2.5). The agreement between the analytic predic- 
tions (summarized in Fig. [2]) and the simulation of the 
EW process in an annealed random network is good in 
the asymptotic large system-size limit. 

Figure [21] shows numerical results (from numerical di- 
agonalization) for both rough phases in quenched net- 
works. For the quenched case, the recurrent /rough phase 
II is somewhat "degenerate" in that it collapses onto 
a single point, a— 2 (Fig. [2]). Here, the asymptotic 
analytic results indicated that the exponent of the di- 
vergence depends continuously on p, G(0) ~ £7(p)-2 
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FIG. 19: Numerically obtained finite-size propagators G(0) 
vs. p for a—lA, q—1 (transient/smooth phase II) for various 
system sizes. Solid symbols represent data obtained by exact 
numerical diagonalization averaged over 100 realizations for 
the quenched network. The solid straight line represents the 
slope obtained from the perturbative analytic calculations in 
the L-^oo limit. Open symbols correspond to the steady- 
state averaged propagator on the annealed network with the 
straight dashed line being the analytic result for the annealed 
system in the L—>oo limit. 

[Eq. ([93]) ]. Although the self-consistent formula pre- 
dicts this feature qualitatively (not shown on the plot 
for the sake of clarity), the actual exponent appears to 
be strongly affected by higher-order corrections. In the 
recurrent/rough phase I, the width diverges as G(0) ~ L, 
just as in a network without long-range links. 

E. Short- and long-time behavior of the expected 
number of returns of a random walker on PL-SW 
networks 

In experimental and numerical setups, it is impor- 
tant to understand the behavior of the observed quan- 
tities for hnite times and finite system sizes in order 
to interpret the results. Because of finite system sizes, 
one can observe different crossovers in the behavior of 
the expected number of returns of a random walker. 
One of the crossover times, T x , can be predicted by 
studying the integral representation of G(0,cj = 1/T) 
in Eq. p62[) : determining whether, at timescale T, the 
behavior of the integral is dictated by the k 2 term or by 
E(fc), the self-energy generated by the random links. The 
other crossover happens at times, T xx , when T xx /L d w 
G(0,w = 1/T XX ). For T > T xx , the scaling of F(T) is 
dominated by finite-size effects. 

As an example, let us have a closer look at tran- 

-2 

sient phase II. Before T x oc p 2 ~ a , the expected num- 
ber of returns scales as that of a regular one-dimensional 
walker; for later times, typically, the walker escapes from 
the vicinity of the origin through some long-range link; 
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FIG. 20: System-size dependence of the width obtained by 
numerically integrating the EW process on annealed PL-SW 
networks [Eq. (|163[) ] (for q=l and p=l) at a representative 
point in the transient/smooth phase II (a=1.4) and in the 
recurrent /rough phase II (a— 2.5). The solid line corresponds 
to the analytically-predicted asymptotic scaling in the recur- 
rent/rough phase II, G(0) ~ L a ~ 2 . The dashed line indicates 
the exact scaling behavior in a network with no long-range 
links in one dimension, G(0) ~ L, as a reference. 



though, after T xx oc Lp^^ , the walker starts to return 
"from the perimeter" of the network due to the finite size 
of the system. 

In Fig. [HJa), we show the sketch of the scaling pre- 
dictions for F(T) as discussed above. In Fig. I^T b). we 
show the disorder- averaged F(T), obtained by employ- 
ing Eq. (|16[) with the eigenvalues from exact numerical 
diagonalization, for each realization of the network, and 
averaging over 100 realizations. The two plots compare 
very well. 



F. Finite-size behavior and scaling functions 

Up to this point, in the analytic forms, we mostly con- 
sidered the infinite system-size behavior of the width, al- 
though it is clear from the plots, Figs. IT51 and [T91 that in 
the limit of weak long-range interactions the divergence 
of the width is cut off by system size effects. This ob- 
servation is plausible since in the case p — ► (or q — > 
obviously) the width should not have a stronger diver- 
gence than the pure (in this case one-dimensional) lat- 
tice. Therefore, one can conclude that for p = (the 
limit of a regular one-dimensional network) [(u> 2 )] ~ L, 
while for p ^ 0, in the infinite system-size limit, it ap- 
proaches a constant, [(w 2 )] ~ l/v£ = £ [see Fig. [23la~)]. 
Thus, the finite-size behavior of the average width can 
be expressed as 
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♦- ♦ a=2.0, p=0.4 
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FIG. 21: Comparison of system-size dependence of the 
width calculated by exact numerical diagonalization (all data 
for q=l) for representative a-s from different phases of the 
quenched network as predicted in Fig. [3] Data in the recur- 
rent/rough phase I (o symbols) indicate the linear scaling 
with L (data is shifted up for clarity). Filled symbols rep- 
resent the scaling behavior in the recurrent/rough phase II 
(a— 2) for various values of p, indicating the p-dependence of 
the exponent. The solid line represents the linear system-size 
dependence to guide the eye. Data in the transient/smooth 
phase II (□ symbols), indicating a finite G(0) in the L^oo 
limit, are also shown for comparison. 



where f(x) is a scaling function such that 



x if x <C 1 
const, if x 1 



(169) 



[(w 2 )} = LfWL) , 



(168) 



For the hard SW network, £ = q + ... hence £ = 
1/y/q- The scaled numerical data, [{w 2 )]/L vs q~ x l 2 jL 
in Fig. l23Tb). shows good collapse, as suggested by 

Eq. Am 

In Fig. [2H for quenched plain SW networks, a simi- 
lar finite-size correction can be observed. From the per- 
turbative calculations Sap 2 but the constant factor 
and higher-order terms were undetermined. In conclu- 
sion, £ oc p -1 . Using this knowledge and doing the same 
scaled plot as for the hard network, a good data collapse 
is observed, though due to our lack of knowledge about 
the accurate behavior of £ the scaling function is not as 
precise as in the hard case. 

Strictly speaking, power-law networks do not have a 
correlation length like plain SW networks (a — 0). Still, 
something can be said about the finite-size behavior, as 
observed in Fig. [19J In the general case, there are two 
competing terms determining the large scale behavior of 
the Green's function, - the diffusion through the regu- 
lar links, A, and through the random ones, A rnd . At 
small length scales, the number of the random links is 
small, their contribution to the ensemble average is negli- 
gible, and the behavior is close to that of a d-dimensional 
unperturbed system. The effect of the random interac- 
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FIG. 22: (a) The analytic prediction for the expected number 
of returns and (b) its numerical verification (fat lines) in tran- 
sient phase II on a quenched PL-SW network with a = 1.4, 
L = 6000, and q = 1. Thin lines represent numerical re- 
sults obtained by excluding the effect of the uniform mode in 
Eq. (|16|l responsible for the "finite-size" behavior for large T 
values. 
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FIG. 23: (a)The system-size dependence of the width for 
different random connection strengths (q-s) on hard SW net- 
works, (b) The data of (a) collapsed to a single curve approx- 
imating E to first order in q and rescaling the axes according 
to Eq. (|168[l . Solid-line segments correspond to the asymp- 
totic small-x behavior of the scaling function f(x) given in 
Eq. (fT69ll . 



tions takes over when the self-energy generated by these 
random links start to dominate the infrared behavior of 
G(k), which is when k 2 , » E(fc x ). Having these different 
regimes and crossovers, one can construct a scaling func- 
tion of G(0) for finite systems similar to the one of plain 
SW networks: 



[(w 2 )}=G(0) = Lf(L x /L) 



(170) 



where L x = k^ 1 oc p 2 -<* in transient/smooth phase II. 
The collapse of the scaled numerical data in Fig. [25] sup- 
ports the validity of the above construction, though slight 
deviation can be observed from the predicted straight line 
for large system sizes. 



G. Disorder-averaged two-point function 



Finally, let us investigate the spatial behavior of the 
disorder-averaged GF, G(l) for quenched plain SW net- 
works. From Eq. ([77]), S(fc) oc p 2 , i.e. the GF is massive. 
In this case, from Eq. (|D3p . 



[G(l)] 



2VS 



(171) 



In Fig. [26] we compared the above analytic result with 
ones from exact numerical diagonalization and found 
good agreement up to finite-size effects at large distances. 
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FIG. 25: The scaling function Eq. (fT70j) for quenched PL-SW 
networks with a = 1.4 obtained by exact numerical diagonal- 
ization of the diffusion matrix. The straight line represents 
the slope predicted by the finite-size scaling argument. 



1 B ^**^r***(f***t* 




-A j> .1 n 1 9 i 

10 10 10 10 10 10 10 10 

ipL)' 1 

FIG. 24: (a) Width as a function of the system size with 
different crossovers for different p-s for quenched plain SW 
networks, (b) Collapse of the same data to yield the scaling 
function Eq. (|169[) . The solid-line segment correspond to the 
asymptotic small-x behavior of the scaling function f(x) given 
in Eq. (TIM)) . 
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FIG. 26: Disorder-averaged two-point function as a function 
of the separation I in simple SW networks for p — 0.10 and 
three system sizes. The solid line corresponds to the expo- 
nential decay given by Eq. (|17ip . The inset shows the same 
data in linear-log plot. 



VI. SUMMARY 

The addition of random long-range links to a regular 
d-dimensional network, producing a SW network, leads 
in many cases to a crossover to mean-field-like behavior 
[HI, effectively becoming equivalent to averaging over 
the long-range links in an annealed fashion. Here, we in- 
vestigated diffusion processes on distance-dependent SW 
networks, where it is not the case (in low dimensions), 
and the contrast between the quenched and annealed sys- 
tems is strong. The results are summarized in the phase 
diagrams of Figs. [2] and [31 with dramatically different be- 



havior in the quenched and annealed case. The results 
for the quenched network were obtained by self-consistent 
perturbation theory. Although the formalism was set up 
for the propagator of the diffusion operator, it can be 
applied to other processes where the mean-field behavior 
is violated or, if the mean- field theory holds, to obtain 
higher-order corrections to it. We also demonstrated the 
validity of the asymptotic theoretical results by providing 
extensive numerical data based on the spectral decompo- 
sition of the coupling matrix and numerical integration 
of the stochastic EW process. 
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APPENDIX A: THE FOURIER TRANSFORM OF 
THE DIFFUSION OPERATOR 



In a general case, the real space behavior of the aver- 
aged diffusion operator (or self-energy) is 



[A 



Jr,r' 



p if r = r' 

AT |r-r'| c 



=p^— ifr^r' 



(Al) 



Multiplying it with a Fourier mode, e zkr , 

V[A™?]e lkr ' = p^e lkr -^V e ,;kr ' 1 l — 

r' r'^r 1 

_ P_ ikr \p 1 (-I _ p*(r'-rh 

r'fr 



P_ ikr ST^ 1 / 1 „ikr"N 

JV ^ |r"l" 1 

r"^0 



(A2) 



The Fourier vectors are eigenvectors of the matrix with 
eigenvalues 



[A™ d ](k) = il 



P (e' kr " - 1) 



(A3) 



r"^0 



or, in the continuum limit, 



p ikr _ i 

[A™ d ](k) = jjjd d r- 



(A4) 



Next, the scaling properties of [A rnd ](k) are calculated, 
when a <C 1/fc <C L. For this reason Af is also approxi- 
mated by its integral-form limit 



L d - a if a < d 



^« / ^"^i^ HL/a) if a = d 
1 ■' " I a d - a if d < a 



(A5) 



The approximations of the formulae below follow those 
of Appendix O Eqs. |C3j and jC9|: 



d d i 



cx k' 



1 — cos(kx) 
l ' k dxx 2 



dx 



L dx 

1/fc * 

L 1 -" ifa<l 
fc"" 1 if 1< a < 3 
fc 2 a 3 "" if 3 < a 



(A6) 



In higher dimensions, 



a zkr 



d d r- 



drf 



d-i 



„ikr cos t 



dn- 



oc ft' 



oc < 



drr 



L d-a 

ln(fcL) 

k a - d 

fc 2 |ln(fca)| 



+ 

'1/fc r 

if a < d 
if a = d 

ifd<a<d + 2 
if a = d + 2 
if d + 2 < a 

(A7) 



Therefore, at leading order, assuming a = 1 without loss 
of generality, 



[A rn<1 ](k) cx ^ 



p if a < d 

p(i + ht) if « = rf 

pfc a ~ d if d < a < d ■ 

pfc 2 |lnfc| ifa = d + 2 

t pfc 2 + ... if d + 2 < a 



(A8) 



In Subsection IIIIFI the behavior of the Fourier trans- 
form of S(r) = — 4 , s „ 7 , ? is needed in two dimensions 

v / n r 4 ln^ {r/aj 

when > /x > 1: 



£(ft) 



s 2tt 



oc — ft 

n 



(Jo(fcr) - 1) 
r 4 ln M (r/a) 
i/fc 1 



rdr 



r hr" (r / a) 



dr 



-hf In ~ M (fca) + 



hr"(fcc 



(A9) 



where ci(fi) is a constant, the result of the Fourier trans- 
formation. 

For the present purposes of this work in most cases 
the scaling behavior of [A rnd ](k) is sufficient to know. 
For more accurate calculations, one has to cope with the 
divergences of the integrals at their limits (for example, 
by introducing counterterms to the integrands) so as to 
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obtain the constants in the formulae. Since it is necessary 
in Subsection IIII El let us calculate E(fc) (or equivalently 
[A r ](fc)) with more scrutiny in one dimension when 2 < 
a < 3: 



E(fc) 



n J a r-y n J a {irp 

1 - e- kr , 



dr 



2 n Sm (r) 
2^fc^ 1 sin(| 7 )r(l- 7 ) 



(A10) 



where 9? is the real part of a complex function. In the 
above steps, we transformed the complex contour integral 
from the real axis to the imaginary one between ia and 
iL, and then moved these limits to and oo since the 
integrand was convergent at these limits and the main 
contribution of the integral comes from fc values oc 1/r 
while a <C r -C L. 



APPENDIX B: THE SCALING PROPERTIES OF 

G(0) AND G(0,u = 1/T) 

Most of the results of this paper are derived for G(r — 
0), or G(0) for short, and G(r — 0, u). Here, their scaling 
properties are obtained for general E(fc)-s when 



G(fc) 



and similarly 



G{k,u) 



k 2 + E(fc) 



1 

fc 2 + E(fc)+w 



(Bl) 



(B2) 



First, let us investigate G(0) 



G(0) 



d d k 



1 



(27r) d fc 2 - 

Sd-i f 1/a 

{^Y Jl/L ' 



■E(fc) 

k d ~ l dk 

: 2 + E(fc) 



(B3) 



Since, in the present work, our interest is only in the 
scaling properties of G(0), the following approximations 
are made: for different k values G{k) is dominated by 
different terms in the denominator of Eq. (|B1[) . as it is 
demonstrated in Fig.[^3for the specific case when E(fc) = 
sk a ~ d . The crossover between the different regimes is 
determined by fc x defined by fc 2 = E(fc x ). From now on 
we assume that 



E(fc) 



0. 



As can be seen in Fig. [23 



k 2 + E(fc) 



(B4) 



(B5) 



if fc x < A;. 



Therefore, 



; fcx k d ^dk f 1/a k d - l dk , N 

(W) ^ /L ^w + L — (B6) 



log(G(k))A 1/L 




log(k) 



log(G(k,w))n 



1/L 
I 



k x 1/a 

H h- 



log(k) 




G(k,w=l/T> 



FIG. 27: A sketch of the approximations used in Eq.-s ()B5[) 
and (|BTT) to obtain the scaling of G(0) and G(Q,u) = 1/T) 
for the specific case when £(fc) = sk a ~ d . 



Next, G(0) will be calculated for specific E(fc) as ap- 
pear in the paper: 
(i) If E(fc) = s, fc x = Vs and 



G(0) oc 



fc d -Mfc /- 1/Q fc d -Mfc 



'1/L 



s" 1 / 2 if d = 1 
= <( ln(sa 2 ) if d = 2 
a 2 - d if d > 3. 



(B7) 



(ii) E(fc) = sk a - d where d < a < d + 2^, k x = s^= 
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and 



G(0) oc 



fex k^dk 
/L sk a ~ d 



l ' a k^dk 



k 2 



-\k 



2d-Qi*-'x 
h/L 



[A' 



-1 

S 3- Q 


if d = 1 and a < 2 


s- 1 ln(sL) 


if d = 1 and a = 2 


l L a-2 


if d = 1 and 2 < a 


\n\sa a - 4 ) 


if d = 2 


a 2 - d 


if d> 3. 



(B8) 



large. This tells us that the 1/T term is irrelevant in 
the integrand while calculating the scaling properties of 
G(0,w = 1/T), the low k behavior is not affected by 1/T 
(for calculational purposes one may interpret it as kx — 
1/T). Another case one must be careful with, is when 
At > k x : in this case E(A) becomes irrelevant in the 
determination of the integral since this case corresponds 
to the scenario that 1/T cuts off the low k behavior even 
before £(A) would have contributed to the integral at all 
(see Fig. |27| . Keeping these in mind the calculation of 
the scaling of G(0, uj) is straightforward. 

APPENDIX C: THE SCALING PROPERTIES OF 

(G(0) - G(r)) 



(hi) E(A) = sk 2 \ ln(Aa)|, A x = ie 1 / 8 and 

rfcx k d - l dk f 1/a k d ' x dk 



G(0) oc 



l/L sA 2 |ln(Aa)| J k 



k 2 



1 



[(d- 



2)-'k d -*\Hka)\-^ L 



if d = 1 



[A 



d-2il/a 

fex 



s In L 

MnlnT if d = 2 

a 2 ~ d if d > 3. 

(B9) 

(iv) £(fc) = sA 2 , £(A) is negligible to the A 2 term in the 
integral and the calculation is straightforward 



L if d = 1 
G(0) oc ^ ln(T/a) if d = 2 



(BIO) 



a 2 - d ifd>3. 



Since (G(0) — G(r)) appears frequently in the calcu- 
lations, its long-distance properties are calculated and 
summarized in this section. 



(G(O)-G(r)) = J 



1 - e ikp d d k 
k 2 + E(A) (2n) d 

dn 



1/a dkk d - 1 f 1- e zkrcosd 



l/L 



(2ir) d J A 2 + £(A) 



(CI) 



In one dimension, 



In order to extract the long-distance behavior of this in- 
tegral, the approximation 



Note that a similar approximation can be done as 
above for G(0,u = 1/T) (see Fig. |2?]). Introducing A T , 
defined as 1/T = E(A T ), 



T if k < k T 

if A T < k < fc x (Bll) 
if A x <C A. 



A 2 + E(A) + 1/T 



Though, one must be careful because, for some £(A), the 
definition of At cannot be satisfied for any A-s (while 
l/L < k < 1/a): for example, if E(A) = s and T is 



(1 — cos(z)) 



z 2 if z « 1 
const, if z >> 1 



(C3) 



is used . 

In higher dimensions, d > 1, 



(G(0) - G(r)) 



I /■ i „ikr cos 9 

■ dk dn \„ a^ 1 



From Appendix ID] 



A 2 + E(A) 



(C4) 



' sm d - 2 8dSl = J^iz)^ 2 * ir 1 ' 2 r 



d- 1 



(C5) 
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The asymptotic behavior of J v (z) is 
and 



Therefore, 



T(v+l) \2J \ L 2r(^+2) \2> 

'^-cos (z- f - *£) 

■kz V 4 2 / 



j o 



d0 



z 2 ,^ r (ff) 



r 



if j» « 1 
if z » 1 



(C6) 



(C7) 



if z << 1 



if z » 1 



(C8) 



Since, in the present work, we are only interested in the (ii) = sk a ~ d where d < a < d + 2 [9J], fc> 
scaling properties of the quantities investigated, it is suf- and 
ficient to conclude that after performing the solid angle 
integrations 



= s 2 + d - 



{e ikrcosd -1) dnoci (fcr)2 iffcr<<1 (C9) 
const, if 1 <C kr. 



Next, using the above approximations, let us calculate 
the scaling of (G(0) — G(r)) when r is large, for spe- 
cific forms of appearing in the calculations. The 
crossover wavenumber k x is defined by fc 2 — S(fc x ). 
(i) S(fc) = s, fc x = \fs and 



(G(0) - G(r)) oc r 2 —k^dk 
hi 



'l/L 



k^dk f 1/a k^dk 



= < 



l/r s Jfcx 

S-V2 ifd=l 

ln(sa 2 ) if d = 2 

a 2 - d if d > 3. 



A: 2 



(CIO) 



(G(O)-G(r)) cx r 2 / ' ' -J— ' -//, 



l/L 



k^dk f 1/a k^dk 



'l/r 
/•" — [fc 



sfc c 



1 r - 2+2d-aql/ r i 1 r;„2d-ai fc x 



[fc 



S 

d-2]l/a 

fev 



J l/L 



fc 2 

k 2 U , 



= < 



s»-° if d = 1 and a < 2 
i ln(rs) if d = 1 and a = 2 
\r a - 2 if d = 1 and 2 < a 
ln(sa Q - 4 ) if d = 2 
if d > 3. 

(Cll) 



,2-d 
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(iii) S(fc) = sfc 2 |ln(fca)|, fc x = \e x l s and 

(G(O)-G(r)) cx r 2 / 7 _ r7rr ^l_fc<'- J ,//,- 



l/L sfc 2 |ln(fca)| 



k d ~ l dk 



1/r sfc 2 |ln(fca)| 7 fex fe 



1/Q fc^dfc 



„2 x r;„d-2 



[r- 2 |ln(fcc 



-liiA- 

J l/L 



+ - [(d-2)- 1 fc d - 2 |ln(fca)r 1 ]i 



/r 



[A- 



d-2il/a 

1 1 



if d= 1 



5 r ln(r/a) 

oc <( ±ln(sln(r/a)) if d = 2 (C12) 
a 2 - d if d > 3. 



(iv) = sfc 2 ln^ ^ (ka) in two dimensions while < 
' and 

1/Q 2tt(1 - Jo(fcr)) kdk 



fx < 1. fc x = ^ esl/<A1 * an( i 



G(0) G(r) " ' /L Jfe2 + afc2ln 1 -''(jfeo)(27r) s 
1 /■ fcx dk 

OC - / ; h ... 

sj 1/r fcln Wl (fca) 

(C13) 

where C2(/x) is a constant, the result of the inverse- Fourier 
transformation. 

(v) £(fc) = sfc 2 , is negligible to the k 2 term in the 
integral and the calculation is straightforward 

r if d = 1 
(G(0) - G(r)) oc <( In(r/o) if d = 2 (C14) 
a 2 - d if d > 3. 



APPENDIX D: G(r) FOR MASSIVE MODELS 

In the literature massive GF are usually appear in their 
/c-space representation: 



G(k) 



1 



k 2 + m 2 



(Dl) 



In such models the correlation length, in terms of the 
mass, is 



£ = m 



(D2) 



Here, the real space behavior of massive GF-s is reviewed. 
In one dimension, 



G(x) = 



dk e ikx _ 1 
2-7T k 2 + m 2 2m 



(D3) 



assuming x > without loss of generality. The integra- 
tion was performed using contour integration by complet- 
ing the integral contour with a semi-circle in the upper 
half complex plane and using the residuum-theorem with 
a pole at ko = im. 

In higher dimensions, d > 1, 



G(r) 



d d k 



ikr 



{2ir) d k 2 + m 2 
dkk^ 1 1 
(2ir) d k 2 +m 2 



(D4) 

dne ikrcose . (D5) 



The integration over the solid angles, dO, can be per- 
formed 



dfle 



ikrcosO n I JO~'~d—2n ikr ( 

= bd-i / aft sin ae 
Ja 



where 



S d = dtt = 



2ir d / 2 
r(d/2) 



(D6) 



(D7) 



the surface of a d-dimensional unit sphere. From one of 
the integral representations of the Bessel functions , 



J v {z) 



7T 1 / 2 T(P+ 1/2) Jo 



dOsin^Oe 



2v c\Az cos i 



(D8) 



and 



t u+i 



As a result, 



G(r) 



■ Jd-2 (kr) 



(D9) 

- 1/2 r(^ + i) 



(27r) d k 2 + m 2 ~ ^ry 

g d _ 1 ^/ 2 r(^i)2^ ^ 



{2n) d ~ X 
I dk— -Jd_ 1 (kr) 



k 2 + m 2 2 ' 



Sd ^/ 2 T( d ^)2^ ^ d 
— r 2 m 2 Kd_ 1 (mr) 



(27T) 



d-2 

m 2 



7* 2 2 

The asymptotic behavior of K v {x) is: 
If i/ = 0, 



(D10) 



If v j= Q, 



-7-ln(|) + ... if « « 1 
/^ e -«(l + 0(I)) if«»l 



r(H) /2\l"l 



if u << 1 



(I) 

£ e -«(l + (I)) if tt »l 



(Dll) 



(D12) 
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In some examples of the main text, the real space 
behavior of the one-dimensional masslcss propagator is 
needed to know. In Fourier space, G°(k) = k~ 2 . The 
real space behavior of it can be obtained by 



G(x) 



— f 



2n/a e ikx dk 



U2 

t/L K 
■y roo ^ikx 

2lT J2n/L 



1 

2^ 



° lkx dk 



l)dk 



2-rr/L 
1 



k 2 



k 2 



dk 



2-k/L 



(D13) 



First, since our interest is in the case a « r « I and the 
integral is convergent at the upper bound so a — > can 
be taken. Second, a counter term is introduced to the 



integrand so as to make the integral convergent at the 
lower bound too. Next, the L — > 00 limit can be taken 



G(x) « — 



1_ f°° (e lkx - l)dk 

1 r 



2tt 7 k 2 

00 (e* K - l)dn 

2^ Jo 



+ L 
+ L 



In conclusion 



G(x) w L — x'f , 
where 7 = J °° ^—^2 ^ K is a constant. 



(D14) 
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[1] B. Kozma, M.B. Hastings, and G. Korniss, Phys. Rev. 

Lett. 92, 108701 [4 pages] (2004) 
[2] B. Kozma, M.B. Hastings, and G. Korniss, Phys. Rev. 

Lett. 95, 018701 [4 pages] (2005). 
[3] R. Albert and A.-L. Barabasi, Rev. Mod. Phys. 74, pp. 

47-97 (2002). 

[4] S.N. Dorogovtsev and J.F.F. Mendes, Adv. Phys. 51, pp. 

1079-1187 (2002). 
[5] M.E.J. Newman, SIAM Review 45, pp. 167-256 (2002). 
[6] S. Boccaletti, V. Latora, Y. Moreno, M. Chavez, and D.- 

U. Hwang, Phys. Rep. 424, 175-308 (2006). 
[7] D.J. Watts and S.H. Strogatz, Nature 393, pp. 440-442 

(1998). 

[8] M.E.J. Newman, J. Stat. Phys. 101, pp. 819-841 (2000). 
[9] S. Milgram, Psychology Today 1, pp. 61-68 (1967). 
[10] P. Erdos and A. Renyi, Publ. Math. Inst. Hung. Acad. 

Set. 5, pp. 17-61 (1960). 
[11] B. Bollobas, Random Graphs (Cambridge University 

Press, Cambridge, 2001) 
[12] M.E.J. Newman and D.J. Watts, Phys. Lett. A 263, pp. 

341-346 (1999). 
[13] R.T. Scalettar, Physica A 170, pp. 282-290 (1991). 
[14] M. Gitterman, J. Phys. A 33, pp. 8373-8381 (2000). 
[15] A. Barrat and M. Weigt, Eur. Phys. J. B 13, pp. 547-560 

(2000). 

[16] M.A. Novotny and S.M. Wheeler, Braz. J. Phys. 34, pp. 

395-400 (2004). 
[17] A. Chatterjee and P. Sen, Phys. Rev E 74 036109 (2006) 
[18] B.J. Kim, H. Hong, P. Holme, G.S. Jeon, P. Minnhagen, 

and M.Y. Choi, Phys. Rev. E 64, 056135 (2001). 
[19] R. Monasson, Eur. Phys. J. B 12, pp. 555-567 (1999). 
[20] S. Jespersen, I.M. Sokolov, and A. Blumen, Phys. Rev. 

E 62, pp. 4405-4408 (2002). 
[21] S. Jespersen and A. Blumen, Phys. Rev. E 62, pp. 6270- 

6274 (2002). 

[22] E. Almaas, R.V. Kulkarni, and D. Stroud, Phys. Rev. 

Lett. 88, 098101 (2002). 
[23] M.B. Hastings, Eur. Phys. J. B 42, pp. 297-301 (2004). 
[24] P.E. Parris and V.M. Kenkre, Phys. Rev. E 72, 056119 

(2005). 

[25] F. Jasch and A. Blumen, Phys. Rev. E 63, 041108 (2001). 
[26] J. Lahtinen, J. Kertesz, and K. Kaski, Phys. Rev. E 64, 
057105 (2001). 



[27] J. Lahtinen, J. Kertesz, and K. Kaski, Physica A 311, 
571 (2002). 

[28] G. Korniss, M.A. Novotny, A.K. Kolakowska, and H. Gu- 
clu, SAC 2002, Proceedings of the 2002 ACM Symposium 
on Applied Computing, pp. 132-138 (2002). 

[29] A. Kolakowska, M.A. Novotny, and G. Korniss, Phys. 
Rev. E 67, 046703 (2003). 

[30] Z. Toroczkai, G. Korniss, M. A. Novotny, and H. Guclu, 
"Scalability, Random Surfaces, and Synchronized Com- 
puting Networks" , in Computational Complexity and Sta- 
tistical Physics, edited by A. Percus, G. Istrate, and C. 
Moore, Santa Fe Institute Studies in the Sciences of Com- 
plexity Series (Oxford University Press, 2006) pp. 249- 
270. 

[31] G. Korniss, M.A. Novotny, Z. Toroczkai, and P.A. 
Rikvold, "Nonequilibrium surface growth and scala- 
bility of parallel algorithms for large asynchronous 
systems" , Computer Simulation Studies in Condensed 
Matter Phystcs XIII, Springer Proceedings in Physics, 
Vol. 86, editors D.P. Landau, S.P. Lewis, and H.-B. 
Schiittlcr ( Springer- Verlag, Berlin, Heidelberg), pp. 183- 
188 (2001). 

[32] H. Guclu and G. Korniss, Phys. Rev. E 69, 065104(R) 
(2004). 

[33] H. Guclu and G. Korniss, M.A. Novotny, and Z. Racz, 

Phys. Rev. E 73, 066115 (2006). 
[34] J. Kleinberg, Nature 406, p. 845 (2000). 
[35] I. M. Sokolov, J. Mai, and A. Blumen, Phys. Rev. Lett. 

79 857 (1997) 

[36] P. Sen, K. Banerjee and T. Biswas, Phys. Rev. E 66 
037102 (2002). 

[37] C. F. Moukarzel and M. Argollo de Menezes, Phys. Rev. 

E 65, 056709 (2002) 
[38] M.B. Hastings, Phys. Rev. Lett. 91, 098701 (2003). 
[39] P. Sen and B. K. Chakrabarti, J. Phys. A 34, 7749 (2001) 
[40] P. G. de Gennes, Scaling Concepts in Polymer Physics 

(Cornell University Press, Ithaca, 1979). 
[41] G. Korniss, M.A. Novotny, H. Guclu, Z. Toroczkai, and 

P.A. Rikvold, Science 299, 677 (2003). 
[42] T. Petermann and P. De Los Rios, Phys. Rev. E 73, 

026114 (2006). 

[43] S.B. Laughlin and T.J. Sejnowski, Science 301, pp. 
1870-1874 (2003). 



3G 



[44] J. A. Davis, V.K. De, and J.D. Meindl, IEEE Trans. Elec. 

Dev. 45, pp. 580-589 (1998). 
[45] S. E. Halford and J. F. Marko Nucleic Acid Research 32, 

pp. 3040-3052 (2004). 
[46] Berg and Clas Blomberg, Biophysical Chemistry 4, pp. 

367-381 (1976) 

[47] D. Brockmann and T. Geisel, Phys. Rev. Lett. 91, 048303 
(2003). 

[48] S. Redner, A Guide to First-Passage Processes (Cam- 
bridge University Press, Cambridge, UK, 2001). 

[49] S. F. Edwards, D. R. Wilkinson, Proc. R. Soc. London 
Ser. A 381, 17 (1982) 

[50] A. Nagurney, J. Cruz, J. Dong, and D. Zhang, European 
Journal of Operational Research 26, 120 (2005). 

[51] B.D. Lubachevsky, Bell Labs Tech. J. 5 (April- June), pp. 
134-156 (2000). 

[52] A.G. Greenberg, B.D. Lubachevsky, D.M. Nicol, and P.E. 
Wright, "Efficient Massively parallel simulation of dy- 
namic channel assignment schemes for wireless cellular 
communications", Proc. 8th Workshop on Parallel and 
Distributed Simulation (PADS'94), pp. 187-194 (SCS, 
San Diego, CA, 1994). 

[53] B.D. Lubachevsky, Complex Systems 1, pp. 1099-1123 

(1987) . 

[54] B.D. Lubachevsky, J. Comput. Phys. 75, pp. 103-122 

(1988) . 

[55] R.M. Fujimoto, Commun. of the ACM 33, 30 (1990). 

[56] A.K. Kolakowska and M.A. Novotny, in Artificial Intelli- 
gence and Computer Science, edited by Susan Shannon, 
series in Progress in Computer Science Research (Nova 
Science Publishers, Hauppauge, NY, 2005), pp. 151-176. 

[57] G. Korniss, Z. Toroczkai, M.A. Novotny, and P.A. 
Rikvold, Phys. Rev. Lett. 84, pp. 1351-1354 (2000). 

[58] M. Kardar, G. Parisi, Y.-C. Zhang, Phys. Rev. Lett. 56, 
889 (1986) 

[59] A.-L. Barabasi, H. E. Stanley, Fractal Concepts in Sur- 
face Growth, (Cambridge University Press, 1995) 

[60] F. Y. Wu, J. Phys. A: Math. Gen. 37 6653-6673 (2004) 

[61] E. Lopez, S.V. Buldyrev, S. Havlin, and H.E. Stanley, 
Phys. Rev. Lett. 94, 248701 (2005). 

[62] Z. Wu, E. Lopez, S.V. Buldyrev, L.A. Braunstein, S. 
Havlin, and H.E. Stanley, Phys. Rev. E 71, 045101(R) 

(2005) . 

[63] J.S. Andrade, Jr., H.J. Herrmann, R.F.S. Andrade, and 
L.R. da Silva, Phys. Rev. Lett. 94, 018702 (2005). 

[64] G. Korniss, M.B. Hastings, K.E. Bassler, M.J. Berryman, 
B. Kozma, D. Abbott, Physics Letters A 350, 324-330 

(2006) . 

[65] W. J. Tzeng and F. Y. Wu, |arXiv:math-phys/0601048 

[66] G. Korniss, Phys. Rev. E 75, 051121 (2007) 

[67] P.G. Doyle and J.L. Snell, Random Walks and Electric 
Networks, Carus Mathematical Monograph Series Vol. 22 
(The Mathematical Association of America, Washington, 
DC, 1984), pp. 83-149; |arXiv:math.PR/0001057| 

[68] L. Lovasz, Random Walks on Graphs: A Survey in Com- 
binatorics, Paul Erdos is Eighty Vol. 2, edited by D. 
Miklos, V.T. Sos, and T. Szonyi (Janos Bolyai Mathe- 
matical Society, Budapest, 1996), pp. 353-398; 

[69] A.K. Chandra, P. Raghavan, W.L. Ruzzo, and R. 
Smolensky, in Proceedings of the 21st Annual ACM Sym- 
posium on the Theory of Computing (ACM Press, New 
York, 1989), pp.574-586. 

[70] P. Tetali, J. Theor. Prob. 4 101 (1991). 

[71] L. C. Freeman, Sociometry 40, 35-41 (1977) 



[72] M.E.J. Newman and M. Girvan, Phys. Rev. £69, 026113 
(2004). 

[73] M.E.J. Newman, Social Networks 27, 39-54 (2005). 
[74] R. Guimera, A. Dlaz-Guilera, F. Vega-Redondo, A. 

Cabrales, and A. Arenas, Phys. Rev. Lett. 89, 248701 

(2004). 

[75] B. Danila, Y. Yu, S. Earl, J.A. Marsh, Z. Toroczkai, and 

K.E. Bassler, Phys. Rev. E 74, 046114 (2006). 
[76] L. Harnau, R. G. Winkler, P. Reineker, J. Chem. Phys. 

104, 6355 (1996) 
[77] J. Rammer, Rev. Mod. Phys. 63 781-817 (1991) 
[78] N. Goldenfeld, Lectures on Phase Transitions and the 

Renormalization Group, (Addison- Wesley Publishing 

Group, 1992) 

[79] W. H. Press et al., Numerical Recipes in C (Cambridge 
University Press, Cambridge, 1995), 2nd ed., Sees. 11.2 
and 11.3 

[80] I. S. Gradstein, I. M. Ryzhik, Table of Integrals, Series, 
and Products, (Academic Press, 2000) 

[81] One can consider this Hamiltonian as the energy of a 
given configuration of the surface and the PDF as a fi- 
nite temperature distribution function of the possible sur- 
faces, in analogy with the PDF of a canonical ensemble 
in statistical physics. 

[82] In the literature, this construction is also referred to as 
the pseudo-inverse of the network Laplacian. 

[83] To obtain the above result, one exploits that the 
GF is defined in the space orthogonal to the zero 
mode *°, hence, £. G tj = £. E^" 1 = 

[84] See Eqs. JTJ and (g7J 

[85] These are the diagrams where all the crosses are con- 
nected together by dashed lines. Diagrams with crosses 
connected to multiple sets by dashed lines are not single- 
link diagrams, like those of Fig. [9] 

[86] This procedure can also be considered as an iteration 
process since, if the initial ansatz was incorrect, one can 
take the result of the SCF as a new ansatz of E and apply 
the same steps over and over again until E converges to 
a function that does not change in the iteration process. 

[87] Note that the other two processes which have three GF-s 
between two endpoints cancel each other 

[88] Note that the fact that the MF approximation is self- 
consistent grants only that the universal behavior of the 
quenched and the annealed system will be the same (e.g. 
the scaling with respect to p) though other features, like 
constants in the formulae, can exhibit non-universal be- 
havior. 

[89] This property can be obtained by noting that the two 
constituents of the Hamiltonian, H° and qV, must have 
the same dimension. Also, we know that H° ~ k 2 and 
qV ~ q/L in Fourier space. Since k has dimension 
(length)' 1 , the dimensionality of q is (length) d ~ 2 . 

[90] This is the same criterion as that obtained in Eq. (|140|) . 

[91] The scenario is drastically different for quenched strongly 
heterogeneous (scale-free) networks, where the largest 
eigenvalue of the network Laplacian is of the order of 
the largest degree in the network, which scales with a 
power of the system size N. Thus, the required At be- 
comes progressively smaller for larger system sizes. 

[92] The exponent is expressed in this way because of the for- 
mulas appearing in the calculations of the perturbation 
expansion. Note that < a — d < 2. 



37 



[93] Note that, strictly speaking, the large-z approximation 
is not really accurate: the function oscillates between 
and 1 but, since it is multiplied with a slowly varying 
function in Eqn. (|74[) , the approximation is adequate to 
calculate scaling properties. 



[94] The exponent is expressed in this way because of the for- 
mulas appearing in the calculations of the perturbation 
expansion. Note that < a — d < 2. 



